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NUMERICAL APPROXIMATION 
The book gives a careful discussion of selected methods 
for such numerical operations as tabulation, interpolation, 
differentiation, integration (quadrature), and curve fitting. 
Each process is illustrated by an actual numerical example 
complete with working detail and with notes on features 
of importance. 


The selected numerical processes are embedded in a 
discussion of the approximate representation of tables of 
numbers and of ‘complicated’ functions over limited 
intervals using polynomials of low degree. Such repre- 
sentations form a natural basis for the development by 
finite-difference and least-square techniques of the various 
numerical processes which are introduced. 


An approach of this kind, in which particular numerical 
processes are set against a proper mathematical back- 
ground, provides students with a reasonable perspective of 
numerical analysis as a branch of mathematics and not a 
mere catalogue of formulae. 


‘Anyone who is concerned with numerical computation 
or who is interested in the application of mathematics in 
such a way as to extract from a table of figures the most 
precise information that it is capable of giving will find 
this book useful.’ — Engineer 
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Preface 


NUMERICAL computation plays an important and growing 
role in the physical sciences, engineering and technology, 
and it serves across the whole range from the development 
of new theories to the analysis of experimental data. As 
problems increase in complexity and as high-speed com- 
puting πο ἀνὰ spread so that their use is the rule rather 
than the exception, it becomes increasingly urgent that 
most students in these subjects should have some introduc- 
tion to numerical techniques. 

The effective use of computing machines involves two 
stages. First we need an understanding of the kinds of way 
in which numerical processes can be constructed and put to 
use; and then we must learn the system of programming 
for the particular machine at hand. Programming systems 
are like languages which enable us to describe problems in 
terms acceptable to the machine; at best they are remark- 
ably easy to learn, but in the long run their effective use 
must depend on a proper preparation of problems and 
hence on a thorough understanding of numerical analysis. 

In this short book (and a forthcoming companion 
volume Numerical Solution of Equations) 1 have chosen to 
concentrate my attention on the mathematical aspects, be- 
cause my experience has shown that in the end students 
will master the use of fast machines more quickly if the 
early stages of teaching are separated, so that they are intro- 
duced first to the ideas underlying numerical processes, and 
then to their adaptation to machine use. 

Like any other branch of mathematics, numerical 
analysis must be understood if it is to be of any real value. 
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PREFACE 


For this reason I have discussed relatively few of the many 
formulae available (as long lists of formulae add nothing 
to understanding), and I have avoided the use of opera~ 
tional methods (as these often confuse inexperienced 
students by obscuring the significance of steps in the deriva- 
tion of a formula). 

The examples and exercises are particularly important, 
as it is difficult fully to grasp numerical processes in the 
abstract; they form a real part of the book, and the reader 
is recommended to follow all numerical steps using such 
calculating aids as may be available. Nothing is simpler 
than to construct fresh exercises from any set of tables. 

B. R. MORTON 
The University, 
Manchester 


CHAPTER ONE 


Introduction to Numerical Analysis 


1.1. The need for numerical methods of analysis 

Progress in the physical, engineering and other sciences 
is based on two parallel but strongly coupled streams. On 
the one hand experimental, involving the design of experi- 
ments and the collection of physical data; and on the other 
hand theoretical, through the construction of model theories 
and their use in the interpretation of experimental results 
and in planning new experiments. The results of experi- 
ments are obtained as numbers or as tables of numerical 
values, and most theoretical solutions have also to be ex- 
pressed numerically before they can be assessed critically 
by comparison with experiment. , 

The most important mathematical techniques used in the 
development of model theories are the normal processes of 
analysis, many of which are introduced in other books of 
this series. These are powerful processes which are of the 
greatest possible importance to the applied mathematician, 
but they often produce formal solutions rather than the 
final numerical solutions that are needed for immediate 
display of the physical significance of a solution and for 
direct experimental comparison. For example, it is easy to 
show that the solution to the differential equation 


dy|dx+2xy=1 with y=0 at x=0 18 ya | “eras, 


but this formal solution is of little numerical value as it is 

easier to solve the original equation by a direct numerical pro- 

cess than to carry out a numerical integration of the rapidly 

growing integrand e. (The ‘formal’ solution is numerically 
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useful given tables of e~** and | jenae, but these will have 


been calculated numerically and the use of such tables is 
itself a numerical process; moreover, this particular in- 
tegral is quite likely to have been tabulated through a direct 
numerical solution of the differential equation!) There are 
other types of problem, such as the solution of non-linear 
differential equations, which are of central importance in 
current work on physical theories but which cannot yet be 
handled analytically to any general extent, and numerical 
methods of solution have an essential role to play here if 
theory is not to lag seriously behind experiment. Again, the 
normal methods of analysis cannot be applied directly to 
operations on tables of experimental values or other tables, 
and such operations must obviously be carried out using 
numerical processes. 

Numerical work will often enter as a normal part of the 
mathematical treatment of a problem, and numerical 
methods form an essential part of the standard equipment 
of any applied mathematician. 


1.2. Approximations to functions 

A basic analytic method is the expansion of a ‘compli- 
cated’ function as a convergent infinite series of terms with 
simpler (or otherwise more appropriate) functional form. 
‘These expansions may be regarded as infinite approxima- 
tion processes in which the error may be made arbitrarily 
small by taking progressively more terms into account. 

The most important case is that of Taylor expansion in 
an infinite power series, where the function f(x) is repre- 
sented for some neighbourhood of a chosen point x» by the 
infinite series, ; 


fla) = Fla) + (1 — sf (0) + 5 (0 — 0) (0) + 


... Ὁ A(x - x0)"f)( a9) +... 
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in powers of the displacement «—xo, where f' =4f dx, ..., 
{Ὁ =dnfldx", and f()(xo) is to be calculated at x= xo. Such 
an expansion will be valid within some range of con- 
vergence a(xo)<« <b(xo), and within this range the trun- 
cation error produced by discarding all terms after the mth 
can be made smaller in magnitude than any given positive 
constant by judicious choice of π. The difficulty from a 
numerical point of view is that this convergence may be so 
weak as to be useless in practice, and also it describes the 
ultimate behaviour of the terms in the series whereas 
numerically we want to concern ourselves with only the 
leading terms. Thus, the expansion 


1 1 1 
tan71 w= x — BH +32? — ax! + eee 


is convergent for —1<«*<1; numerically it is very effec- 
tive for small values of x and only two terms are needed at 
x=0-1 for five decimal accuracy, but the number of terms 
needed increases rapidly with «, to seven terms at x=0-5 
for five decimal accuracy and fifty-one terms at x=1 for 
two decimal accuracy. ee | 

The Taylor series can also be written in the terminated 
form 


fla) =flte)+ (2 οὐ + on 2} ταῦ οι + Rey 
with remainder Ry = αττιο- xo)"+1f+1)( £) for some ξ in 


the interval (xo, x). The remainder R, may be regarded as 
the error involved in the approximate representation 


1 
f(x) =f (x0) + (ὦ — %0)f' (#0) Ἐ ....Ὁ τί — ΚΟ)" (20), 
and it may be seen that this is a simple polynomial repre- 


sentation of degree n for the function f(x). But can we 
reasonably expect that such polynomial representations can 
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in practice be used without substantial error caused by 

neglect of the remainder Ry? This question is answered 

partially by Weierstrass’s theorem on approximation by 
polynomials. 

Weierstrass’s theorem* states that a function continuous 
tn a closed interval can be approximated over the whole 
interval by a polynomial of suitable degree which differs from 
the function by less than any given positive quantity ε at 
every point of the interval. As the allowable error of approxi- 
mation ¢ is decreased, the degree of the approximating 
polynomial must be increased. This theorem is of the 
greatest importance because it implies that in practice, 
where we work always to a limited and specified accuracy 
only, we can replace any continuous function over a given 
range by an approximating polynomial without significant 
loss of accuracy. Hence we are justified in assuming a poly- 
nomial approximation for a given set of experimental or 
other values, and also in matching given functions by poly- 
nomial approximations over closed intervals. There re- 
mains one major drawback, in that Weierstrass’s theorem 
gives no indication whatever of the degree of polynomial 
required for a given accuracy of approximation. This is 
typical of the difference in attitude in pure and in numerical 
analysis; the pure analyst is concerned first with the 
general principles of approximation, but the numerical 
analyst wants to know when he can reasonably use quad- 
ratic or cubic approximation! We shall see that in practice 
there are ways of judging the degree of a suitable poly- 
nomial representation even when no direct estimate for the 
remainder R» can be made. 

_ The way is now open for a major change in emphasis, 
We shall adopt unashamedly the idea of simple polynomiat 
approximation to the functions or tabular values involved 
in calculations as a basis for developing numerical processes. 
Moreover, we shall reinterpret the significance of the re- 


* For a proof and various extensions see e.g. Jeffreys and Jeffreys 
Methods of Mathematical Physics, Cambridge. τῶν ᾿ 
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mainder or error term of the terminated Taylor expansion 
as providing a guide to the working range over which a 
polynomial representation of specified degree can be used to 
the accuracy required. We shall thus construct numerical 
processes in which each function is represented by a poly- 
nomial of reasonably low degree (so that the ‘unit opera- 
tions’ may be easily carried out) and in which the range of 
calculation is subdivided into working intervals which are 
short enough locally to preserve the desired accuracy over 
the whole calculation. 

There are also many types of infinite series expansion in 
terms of complete sets of orthogonal functions. The most 
familiar of these are Fourier series* in terms of cosines, 
sines, or a mixture of both, and in each case valid within a 
range of orthogonality of the set of functions involved. 
Other sets of orthogonal functions which may be used for 
series expansion in an analogous manner include Legendre 
polynomials, Chebychev polynomials, Bessel functions, and 
so on. In principle, each of these expansions can be used in 
truncated form as a basis for the development of numerical 
processes in very much the same way as direct or ‘simple’ 
polynomial representation. In practice, some of the most 
useful methods have been developed from truncated 
Fourier series and series of Chebychev polynomials. 
Fourier series methods form the basis of harmonic analysis 
and are particularly useful in handling functions or data 
with natural periodicities. Chebychev polynomial methods 
provide the most rapid convergence of all approximations 
based on special or simple polynomials, and give a special 
degree of uniformity in the distribution of errors of ap- 
proximation over the range of application (in contrast with 
simple polynomials where the error of approximation tends 
to be smailer in the centre and to increase fairly sharply 
towards the ends of the range). It may be noted that a 
similar form of Weierstrass’s theorem holds for Fourier and 
Chebychev representations. 

See Sneddon, — Series, in this series 
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All the methods of approximation which have been 
described can be collected into the single expression 


H(%) = ao go(x) + a191(x) + ... +andn(x), 


and any such finite approximation may be used as a basis 
for the development of numerical processes (though only 
with complete sets of functions can we be sure that the error 
will be reduced progressively as we include successively 
more terms of the representation). This representation is an 
approximate form with error Rp() of the exact terminated 
expansion 


S(%) = aogo(*) + αἱ φι(χ) + ... +andn(*) + Rn(x). 
In order to emphasize the role of approximation in the 
processes we shall develop, the symbol = will be used to 
denote any of these approximate representations, such as 
polynomial approximation with ¢;(x)=x* or Fourier 


approximation with ¢;(x)= ae akx. Many authors use the 


symbol = for all purposes, but this can prove confusing to 
the beginner. In each case the maximum value of R»(x) (if 
it can be estimated) in a range of application provides a 
measure of the accuracy obtainable with that approxima- 
tion, 

A family of numerical processes will be developed in the 
following chapters from the simple basis of polynomial 
approximation, using the general form 


I(x) = ao αικ +a0x2 + 2... + anx”, - 


(with appropriately chosen degree m) for known and un- 
known functions and for sections of tabulated values. 
Numerical processes based on polynomial representation 
are relatively easy to derive, simple to apply, and very satis- 
factory for displaying the principles of numerical approxi- 
mation ; moreover, they require only a very modest mathe- 
matical background, and so form an ideal introduction to 
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the ideas underlying numerical analysis. In contrast, the 
development of methods based on approximation by special 
functions calls for a wider knowledge of analysis; hence in 
spite of their greater current importance, methods based 
on approximation by orthogonal functions will be post- 
poned to a later book in the series. 


1.3. Numbers, continuous functions, tables of values, 
and tabulated experimental results 
These provide the raw material of numerical analysis. 
We shall normally represent numbers in decimal form.* 
Thus the number 
3-14159 


has five decimals (or decimal places, or decimal digits), and 
six significant figures. It is the value of 7, and is to be 
interpreted as being correct to within five units in the sixth 
decimal place; it is impossible to guarantee greater accuracy 
without writing additional decimals, and it is an important 
convention that it should be regarded as dishonest to fill 
decimal places with digits known to be inaccurate (or, 
indeed, with digits of unknown accuracy). This represents an 
approximate value of 7 which is of course an irrational 
number, but irrational numbers cannot as such enter 
numerical analysis since we are forced to approximate to 
them. In the same way we must approximate all rational 
numbers which cannot be expressed exactly in the number 
of decimals with which we choose to work, and for obvious 
ip reasons this number cannot normally be very 
arge. 

Sometimes we have to handle very large or very small 
numbers, such as Avogadro’s number 6-023 x 1023, or the 
charge of the electron in e.s.u. 4-802 x 10-2; and in such 
cases we usually separate the significant figures (which we 
write in the form of a workable decimal number of order 


* This is by no means the only possibility, though it is most con- 
venient for our purposes; modern computing machines use binary form. 
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unity) from the powers of ten involved, and we carry the 
two parts of the number separately through the calculation. 
This may be called floating point decimal form, and may also 
be written 6-023, 23 and 4-802, —10 for these examples. 
It may be noted that the value for the charge of the electron 
implies that it is known correctly to about one part in ten 
thousand (to be precise, 5 parts in 48020). 

Continuous functions could in principle be represented 
numerically by curves showing their value at every point, 
but for practical purposes the only way in which most 
functions can be represented accurately and conveniently is 
by a table of values at specified intervals of the indepen- 
dent variable.* And of course we are only too familiar with 
such tables of logarithms, sines, cosines, exponentials and 
so on. In these printed tables each entry represents a true 
value, correct to within half a unit in the last digit unless 
someone has blundered, of a precisely defined function. (We 
must accept the fact that other people make mistakes just 
as we do ourselves, and quite a lot of mistakes appear in 
print!) These functions obviously exist between tabulated 
values, and we shall develop the process of interpolation to 
calculate such intermediate values. However, we cannot 
tabulate the functions more closely and many functions are 
tabulated scarcely at all because of the sheer magnitude of 
the task involved. Thus, a table of five-figure logarithms 
for the numbers from 10 by intervals of 0-1 to 100 occupies 
two pages with 900 main entries; a similar table of a func- 
tion of two variables would take two large books of 900 
pages each; one of three variables a room of books; one of 
four variables (such as the hypergeometric function) a 
library of more than one and a half million volumes! The 
evaluation of many functions can be carried out very 
quickly on fast computing machines. | 

Tables of experimental results differ from tables of functions 

* But see Jahnke and Emde, Tables of Functions, Dover for some 


beautifully drawn curves which are extremely effective in displaying the 
broad types of behaviour of rather  ebimat functions. 
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in that each entry is likely to include experimental error, and 
although we can place probable bounds on the magnitude of 
the error we can do nothing to improve the accuracy of indi- 
vidual results without repeating and refining the experi- 
mental measurements. There are two kinds of experimental 
error: systematic and random. Systematic errors may be due 
to a persistent fault in the performance or faulty design of 
the experiment; they can seldom be detected by direct 
analysis of the experimental results as each measurement 
suffers, and they are less a part of the analysis of the results 
than of experimental design. Random errors are often ap- 
parent on inspection, specially when the results are pre- 
sented graphically and show ‘scatter’; they can be kept at a 
low level by careful work but they cannot be eliminated, 
and they set an obvious limit to the accuracy which can be 
obtained in calculations using the experimental results. 

A common way of approximating a set of experimental 
results is to plot them on suitable logarithmic scale or other 
graph paper in an attempt to produce a linear relationship ; 
a ‘line of best fit’ can then be estimated by eye, and will 
serve to define a functional approximation to the original 
results. When the scatter is appreciable, there is some temp- 
tation to use one of the smoothing techniques so that the 
adjusted points lie more closely on a smooth curve; how- 
ever, such devices should be used with great caution, as the 
raw experimental results provide our only direct link with 
the physical world and if they do not form an obvious 
pattern this may be because we have taken the wrong 
measurements. We shall mention more sophisticated 
methods of curve fitting in chapter 5, and it should be noted 
that there are statistical methods for assessing the signifi- 
cance of variations in experimental measurements. 


1.4, Errors and accuracy in numerical analysis 
The only useful numerical solution is a correct one— 
though we must be careful how we interpret the word 
correct’, We can usefully distinguish two causes producing 
Β 9 
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wrong numerical answers: errors which arise from the finite 
character of numerical analysis and are inherent in our 
processes, and blunders caused by our own carelessness of 
the failure of machine aids. 

Inevitably we must work with a limited number of sig- 
nificant figures, because our machines have limited capacity 
and indeed we have even more limited capacity ourselves 
in practice! Few numbers will be small integers, and there 
will usually be a rounding off error produced by the rejection 
of all digits after a selected decimal place. Our procedure 
for rounding off will be to preserve unchanged the last 
digit retained on discarding fewer than five units in the next 
place, to increase it by one on discarding more than five 
units, and to round off to the nearest even number on dis- 
carding exactly five units. Some authors always round five 
units in one direction (usually up), but our procedure 
preserves a randomness and so increases the chance of 
mutual cancelling of rounding errors during arithmetical 
operations. Thus, on rounding to three decimals, 


0:5772 +0°577, 3-14159 + 3-142, 
0:0625 + 0-062,  0-1875 -Ὁ 0-188. 


This procedure normally ensures that the rounding error 
will not exceed half a unit in the last digit retained. 
Rounding errors can accumulate in arithmetical opera- 
tions, and it may be necessary to carry an extra decimal 
digit, otherwise unnecessary, rounding off at the end of a 
chain of operations. The probable error will be consider- 
ably less than the maximum error, and can be estimated 
statistically. For example, the maximum error in the sum 
of n rounded-off values is -—-47 units in the last digit and 
occurs only when each value has the greatest possible 
rounding error and all reinforce; in practice the distribu- 
tion of rounding errors is likely to be almost random and 
for large n the probable error in the sum is likely to be less 


than +4Vn/3, with error about +V n/3 in roughly one 
entry in twenty. 
10 
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Frequently we obtain finite numerical processes by suit 
ably truncating an appropriate infinite series. This neglect 
of higher terms introduces a truncation error, which can be 
kept small by retaining a sufficient number of terms 
(though with rapidly increasing labour as the number of 
terms increases), ‘Truncation errors will be present in most 
processes other than those of arithmetic. 

Blunders are lapses from accurate computation due to 
carelessness of the computer or failure of the machine. In 
principle they are avoidable, but in practice they are almost 
impossible to eliminate; hence we must be constantly on 
guard, and the checking and rechecking of calculations 
must become a part of our habitual practice as computers. 
For example, in adding a column of figures we should 
always do it twice, once from the bottom and once from the 
top; if the sums do not agree we must repeat until two suc- 
cessive checks do give the same answer. Whenever possible, 
checks should be made using different methods of calcula- 
tion, because if we have made a mistake once we are quite 
likely to make it again in a repeat calculation; careful com- 
putational design will often produce formulae which are 
not used directly in the calculation but must nevertheless 
be satisfied, and these provide valuable checks. We should 
eee finish a problem with two consecutive checks that 

gree. 

Beginners often try to laugh off this emphasis on blun- 
ders, but they learn how easy it is: to transpose digits, 
taking 7239 in place of 7329; to misread repeated digits, 
taking 72239 in place of 72339; to misread tables; and to 
overlook signs, especially near sign changes. They find that 
after identifying a mistake they will make a wrong correc- 
tion. These and other blunders can be kept to a minimum 
only when the computer learns to work with an easy 
rhythm, to write his figures clearly and to lay out his com- 
putation in an orderly fashion. And he will discover his 
remaining blunders quickly if he checks the work stage by 
stage, avoiding checks by repetition; finally he will be con- 

11 
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fident of his answer only when he has completed an over- 
all check. 

Inexperienced computers often carry too many—some- 
times far too many—decimals through their calculations 
from lack of confidence, and from a vague idea that they 
may increase the accuracy of their work. We must under- 
stand clearly that although our model theories in science are 
capable of unlimited accuracy mathematically, they are 
based on approximate physical assumptions and are directed 
towards approximate physical measurements. In many ways 
the bridge of mathematics in between is expendable, and 
there is certainly no merit to be gained from unreasonable 
pursuit of accuracy in our model theories. Indeed the very 
striving for excessive accuracy defeats its own purpose; for 
the computer, burdened by excessive length in his calcula- 
tions, will quickly tire and err. We shall do well to recognize 
as the mark of the experienced computer the fact that he 
plans his problem so that at each stage he carries just that 
accuracy required locally to produce a final answer as 
accurate as the experimental measurements and no more 80. 


1.5. Effective numerical computation 

The criterion by which we must judge numerical pro- 
cesses is effectiveness. This involves mathematical considera- 
tions such as accuracy and availability of independent 
checks, but also what may best be described as economic 
considerations such as the time required to obtain specified 
accuracy. When there is a long run of similar calculations 
we shall spend ample time in reducing the problem to the 
best possible form for computation, so that each case may 
be handled quickly and accurately with an over-all saving 
in time and effort. But if we have a single case to compute 
we shall use any already-familiar process, and the calcula- 
tion will be finished in less time than a preparatory reduc- 
tion. We shall, of course, use mathematical tables, slide 
rules, hand machines, and above all else high-speed elec- 
tronic computers whenever ri will be helpful. 

1 


EXERCISES 


A glance at some of the larger texts will show that there 
are very many processes available for each type of calcula- 
tion, such as interpolation or integration, for example. No 
doubt the beginner would wish to see one process of each 
type labelled best, but this cannot be done. There are no 
‘best processes’ ; the most effective processes will always depend 
on the particular problem in hand. There is great scope for 
the use of judgment in numerical analysis in such matters 
as the selection of an appropriate process, the detailed 
design of a calculation, the choice of working interval, 
number of working digits, and so on. Indeed it is this de- 
mand on the capacity for skilled judgment by the computer 
that provides the challenge and much of the interest of 
numerical work. 

In this short book and a companion volume on Numerical 
Solution of Equations we shall discuss the kinds of way in 
which numerical processes can be constructed and put to 
use. We shall introduce few out of the vast array of avail- 
able formulae, but as the beginner enlarges his experience 
he will learn how to select the most suitable from among 
many processes that might be used on a given problem. 


EXERCISES ON CHAPTER ONE 


1. State the number of decimal places and of significant figures in 
each of the following numbers, and (approximately) the accuracy to 
which they are presented: (i) 980-665, (ii) 9-11 x 10-26. (iii) 2.9978 x 10, 
(iv) 1-20 x 107, 

2. Compare the result of rounding off 1:37455 digit by digit until 
only two decimals remain, and rounding off to two decimals in a single 
*peration. 

3. Evaluate and write in appropriate form, giving an estimate of the 
Maximum error in each answer: 

(i) 1-73 —2-16 4+-0-08 + 1-00 —2-23 —0-97 +-3-02, 
(ii) 9°37 x8-91, 
(iii) 7-137/0-279, 

(iv) 29132 +15-164—0-31725. 


4. Compare the sum and a error when the following 
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terms are added and then rounded to six decimals with those obtained 
when the terms are first rounded off and then added: 


“5 3 4+--1000000 +--0083333 + -0011905 - ΟΟΟ2226 + 
ones 1} . + 0000493 -+--0000122 +--0000033. 


CHAPTER TWO 


5. Obtain linear, quadratic and cubic polynomial approximations for 
εἶ near x=0 by truncating the Taylor series appropriately, and use the 
remainder term to estimate (to the nearest 0-1) the range over which each Finite Diffe 

lynomial can be used if the error of approximation is not to exceed Hinite erences 
0-005 (i.e. the approximation is to be correct to two decimals). 


Before we can develop particular processes we must in- 
vestigate the operations that can be carried out directly on 
a table of values. Continuous functions must normally be 
represented numerically by tables of values at specified 
᾿ intervals of the argument, and the limit processes of in- 
| finitesimal analysis can no longer be applied to such tabu- 
lated functions but must be replaced by equivalent finite 
processes. 


2.1. Tables of values 

Most functions in common use are tabulated, and these 
tables of values provide one of our sources of raw material. 
Any scientist should become familiar with selected sets of 
tables: for example Comrie, Chambers’s Shorter Six-figure 
Mathematical Tables is good for ordinary functions, and 
Jahnke & Emde, Tables of Functions is useful for some less 
common ones, To illustrate table notation we may take as 
typical Chambers’s table of ordinary logarithms 1,000 (1) 
10,000, 6D; the table gives in a rectangular array the 
mantissae correct to six decimals (denoted by 6D) for the 
numbers from 1,000 to 10,000 (the endpoints and the tabu- 
lar interval 1 in the argument are denoted by 1,000 (1) 
10,000); decimal points and characteristics are suppressed 
as they convey no additional information of value, and their 
inclusion would make the table harder to work from. Tabu- 
lated values correct to six significant figures are denoted 
by 6S. Another illustration is provided by Chambers’s 
table for δὲ given at values of the argument 0-000 (-001) 
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3-000 (-01) 10-00 (-1) 50-0 to between six and seven signifi- 
cant figures; the entries are in decimal form up to x=15-0, 
and thereafter in floating point decimal form 6D, 7S so 
that the entry for e!5=3,269,017 is printed as 3-269017, 6. 

Some functions change rapidly in isolated regions of the 
argument and slowly elsewhere, and such functions can 
best be tabulated by using different tabular intervals in 
different parts of the range according to the local behaviour 
of the function. The resulting tables are rather more diffi- 
cult to work from than those with uniform interval, and the 
great majority of the tables in common use are based on a 
constant tabular interval in the argument, at least over 
substantial regions (e.g. e* above). A section of such a table, 
say for f(x)=1/x, might be written 


* 30 31 32 33 34 365 36 #37 
1/* :33333 -32258 -31250 -30303 -29412 -28571 -27778 -27027 


This table cannot represent the nominal function 1/x 
exactly; for there will be rounding errors in most entries, 
and indeed the tabulated function is not formally defined 
at all for intermediate values of x. In this case we know very 
well that the tabulated values provide a good approxima- 
tion to 1/x, and that intermediate values can be estimated 
by fitting a polynomial (Weierstrass’s theorem!) and we 
might anticipate (or hope?) that a low degree polynomial 
would be adequate because of the slow variation of 1/x in 
this interval. But what if we are merely presented with a 
table of values? In principle this could represent some patho- 
logical function of pure mathematics, and even if Weier- 
strass’s theorem should be applicable a very high degree 
polynomial might be required for matching; or it might 
represent some perfectly well-behaved periodic function 
with natural wavelength considerably shorter than the 
range of the table. ‘Thus, table 1 gives values for sin 22x 
at three different working intervals h=0-9, 1-0, and 1-05 
for the variable x; it is ‘perfectly clear’ that each of these 
16 
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Table 1 
κ ῦΌ ἃ oh Sh 4h 5h 
09 O —-588 —-951 —-951 —-588 0 
sin2ax<hk=1:0 0 0 0 0 0 0 
1:05 0 ‘309 ‘588 "809 ‘951 1.000 


three sets of values defines a ‘good smooth function’, but 
unfortunately none of these ‘good smooth functions’ bears 
the slightest useful relationship to the original function 
sin 27x! In each of these (rather artificial) cases the work- 
ing interval is so large relative to the periodic length of the 
function that the main periodicity has entirely escaped 
representation in the table. An even worse example is pro- 
vided by the function sin (27/x), for in this case as x de- 
creases, no matter how small a working interval we may 
select it will soon become impossibly large relative to the 
continuously shortening periodic length of the function. 
Clearly there can be no certain way of knowing whether a 
given table provides a proper representation of the (un- 
known) function which it is supposed to represent. The 
interval adopted in printed tables will always be short 
enough to display fully the essential character of the func- 
tion, and it is important that we should always adopt the 
same convention ourselves when giving tables of values; 
and it should be noted that this will usually mean using 
different intervals in different parts of the range according 
to the local behaviour of the function. We shall assume in 
working from tables that the higher order derivatives of 
the function are small enough for the function to behave 
eR ig smoothly’ over the length of the tabular 
int ‘ 


2.2. Finite differences 


_ The simplest operation on a table of values, obtained 
either for a known function or by calculation, is to find the 
difference between each pair of tabular values. The first 
differences are found by subtracting each value from its 
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successor in the table, second differences by repeating a 
similar operation on the first differences, and so on for 
higher orders; these together comprise the finite differences 
of the table. The customary form for tabulation is: 


x fe) = first- second- _ third- 
1/x differences 
80 33333 
— 1075 
31 ‘32258 67 
— 1008 —6 
32 31250 , 61 
— 947 --ὃ 
8: "90808 56 
— 891 —6 
3°4 29412 50 
— 841 -2 
8:8 ‘28571 48 
— 793 -θ 
8:6 ‘27778 42 
— 751 
8:1 ΕἼ ΡΥ] 


Note: (i) that odd-order differences are written between, 
and even-order differences on the tabular lines; (ii) that 
decimal points and leading zeros are omitted in the dif- 
ferences in order to simplify the table; (iii) that we must 
very definitely not omit the signs; (iv) that each difference 
must be multiplied by 10-5 before use. 


Two features stand out in this table. (1) Differences of 
increasing order decrease quite rapidly in magnitude. We 
shall show in §2.4 that this type of behaviour is typical of 
polynomials and of functions or tables which are well-fitted 
for polynomial approximation. (2) The third differences 
show marked irregularity. If, however, we tabulate 1/x to 
six decimals (6D) we find that the retention of the additional 
decimal has produced fairly regular third differences, but 
markedly irregular fourth oO 
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* 1/x first- _ second- _third- _fourth- 
differences 
90 .333838 
| — 10752 
31 ‘322581 671 
— 10081 — 60 
8:2 *312500 611 7 
— 9470 — 63 
3°3 303030 558 3 
— 8912 — 50 
3°4 294118 508 10 
— 8404 — 40 
35 285714 468 0 
— 7936 — 40 
3°6 277778 428 
— 7508 


9.7 “270270 

This behaviour can be seen to arise from rounding errors; 
even though individual rounding errors in f(x) are <4 in 
the last decimal digit it is clear that these errors are likely 
to accumulate, and be increasingly amplified in succes- 
sively higher differences. The greatest error that can be 
obtained is, approximately, 

tabular differences 
error ist 2nd 38rd 4th bth 6th 7th 


—1 
an +2 
+1 —4 
+4 —2 +8 
—1 +4 —16 
-Δ +2 -8 +32 
+1 —4 +16 — 64 
+4 —? +8 — 32 
--1 +4 —16 
-; +2 - 
+1 —4 
+4 —2 
-1 
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These are maximum errors, and the likely errors will be 
smaller as the rounding error will not usually alternate in 
sign; Comrie* has given a working criterion for the ex- 
pected fluctuations which will arise from round-off as 


order of difference ἘΠ 2.8 5 6 
expected limits of error +1 +2 +3 +6 +12 +22, 


and we must learn to tolerate small irregularities of this 
magnitude in our tables of differences. This level of fluctua- 
tions has been called the ‘noise level’ of a difference table. 

The magnitude of the finite differences obviously de- 
pends on the length of the tabular interval. 

Although tables of values play a more important role in 
hand calculation than with high-speed machines, they are 
generally important for the expression of experimental 
results and of the results of a previous calculation. They are 
seldom used for the description of specified functions with 
fast machines, as it is usually much quicker to use a pro- 
gramme sub-routine to evaluate a function where it is 
needed than to feed in a table of numbers. 


2.3. Notation for finite differences 


There are three distinct notations for the single set of 
finite differences of a function f(x) tabulated at constant 
interval h in the independent variable x. If we write 
Xn =xo+nh for the mth tabular point, and f,x=/(xn), then 
we define 

(i) The forward difference at x is ἔμ1-- ἴ written Af, 


hus 
Af(xo) =f(%0 + ἢ) —f(*o) =fi—fo, 
and higher differences 42fy = 4(Afo) = Δῇ -- 4fo =f2—2fi +fo, 
Af. = A(23fs) = =f5— 3. + 38 — fo, etc., can be calculated 
either by differencing lower differences or from the tabu- 
lated values. Note that 


* See introduction to Chamber's Six-figure Mathematical Tables, 
Vol. 1]. 20 


2.3 NOTATION FOR FINITE DIFFERENCES 
A*fi=frst— (i) fore ΕΣ (5) frri-o— ... +(-1)*%f;, 


where (*) are the binomial coefficients. 


* f(x) 
Xo to 
Afo 
ie, fi A*fo 
Af, Af, 
v2 fe Af, 1 
Afr 
xg fs, 
Xn-3 = fn-8 
Pfn-2 
“n-2  fn-2 P*fn-1 
‘i f fn of γι 
An—1 n—-l Γ n 
Vin pee asd 
vn Ii ik = ell εἰ 


(ti) The backward difference at x; is fi—fi-, written 
γί. Thus, for example, 


γ᾿ = 7 (fn —fn-1) =f - 2fn-1 +fn-2. 


The forward and backward differences represent in 
aggregate precisely the same set of numbers (and 
Vft=4ft-1), but forward differences are specially useful 
near the start of a table as they involve only the tabular 
values below the ‘subscript level’ x;, and backward differ- 
ences are useful near the finish of a table where we do not 
have the forward tabular values to obtain all the necessary 
forward differences. 

(iii) The central difference at x; is fi++»—fi—», written 
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Sf. Thus Sfi+%h == +; (+1 —fi = Af; = γί ιν etc. This is the | 


normal notation for use away from the ends of a table, 
where there are plenty of tabular values above and below 
the level χε. 


fi-2 
ofi-), 

{ει δὲῇ,-ἃ 
fia fa 

fi ———> 3 —__> μῃ 
δύ. oy rel 

tin δέῃ δ μ 
ΣΕ 55/443) a 

Fixe 142 
δή i+5/, 

7.3 


Note that 4"/; is the nth forward difference and lies on 
the downwards diagonal from the tabular entry ἤν, while 
the nth backward difference yf; lies on the upwards 
diagonal originating at f;. The nth order central difference 
é"f; lies on the horizontal level χε; thus odd central dif- 
ferences such as δύ lie on half tabular lines and even ones 
like 82f_, on tabular lines. 


2.4. Finite differences of a polynomial 


The finite differences of a general polynomial expression 
can be found only in terms of a specified tabular interval ἢ. 

The finite difference of a constant over any interval is, of 
course, zero. 

The successive forward differences of the power x* are 


Ax* = (0+ h)* —x¥ = khxk-1 + $R(R—1)h2xk2 +... 
+kh*¥-lx + hk, 

Δ" κα A( Ax”) om k(k— 1h? *-8 + k(k —1)(k—2)h3x'-8 
+ vee τ(28 ---2)}}, 


᾿ * , . ῃ ' » 
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Atxk = k(k—1)...(k-i+ 1)htxt-+ + 4ik(k—1) ... 
(k—a)hitgk-tAly 


ΓῚ - ΓῚ ΓῚ 5 * - 


Akxk=k! h*, 
All higher differences 4*+4x*, 1:50 are zero, Thus the ith 
difference of x* is a polynomial of degree k—i for 0<i<f, 
and is zero for i>k. 
The finite differences of an nth degree polynomial 


f(x) = a0 + ax + aor? + 1... +anx” 


follow, and it can be seen that 4‘f(x) is a polynomial of 
lower degree n—12 if 0<i<za, and is zero if i>n. In par- 
ticular, 

A™(ag+ ax+ ... + nx") =ann! h*, 

A®*l(ag+ Qyx+ ... + @nx") τεῦ. 


Thus the nth difference of a polynomial of degree n is a constant 
proportional to h", and higher differences are zero. 


Example 1. Derive finite difference tables for the polynomial 
f(x)=2—3x°+.25 with the two working intervals h=0-5 and 
h - ῦ' 1 ᾿ 

These tables can easily be obtained as: 


% Τα) 4f 4 Δ 41 
0 


2-000 
— 625 
5 1:375 — 750 | 
— 1375 ‘ 750 ; 
1-0 0 
— 1875 750 
15 —1:375 +750 0 
— 625 750 
20  —2-000 +1500 : 
+875 : 
25 «-1:15ὅ ᾿ : 


. ᾿ 
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᾿ f(x) Af a af As 
0 2-000 j hats 
—29 
1 1-971 —54 
— 83 6 
2 1-888 — 48 0 
-- 181 θ 
8 1-757 — 42 0 
—173 6 
4 1-584 —36 
— 209 
5 1:375 


Note: (i) that in each case third differences are constant and 
equal to 3! 18, and fourth differences are zero; (ii) that the dif- 
ferences depend on h and higher differences are systematically 
small only when h is small; (iii) that if tabular values of f(x) 
are rounded off there will be small irregularities in fourth and 
higher differences. 


A comparison of the second table of example 1 with the 
table for 1/x on page 18 suggests that after allowing for 
round-off fluctuation we might very well adopt a cubic 
approximation for 1/x in the range 3-0 <x «3:7. Indeed we 
shall find that whenever the higher differences of a table 
become increasingly small, the function represented by the 
table is well-fitted for polynomial approximation. Large dif- 
ferences may be reduced by reduction of the working in- 
terval h, but in some cases approximation by reasonably 
low-degree polynomials will be possible only for incon- 
veniently small ἢ and other types of approximation may be 
more effective. Thus, it will be found that the values in 
table 1 are ill-suited for polynomial approximation and will 
remain so unless the interval is reduced considerably ; this 
is hardly surprising for a sinusoidal function of wavelength 
small relative to the range of x, and we might expect that 
harmonic approximation would be better. 
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2.5. Detection of mistakes by differencing 

Rounding errors cause only small fluctuations in dif- 
ferences, but blunders will produce much larger upsets. 
Hence a standard and invaluable method for detecting 
mistakes in any table is to form differences and inspect 
those of higher order for irregularity. An error of one unit 
in f(x) (or in a difference), and a pair of adjacent errors 
propagate through the differences as follows: 


0 1 0 I 
0 0 1 0 1 

0 0 1 —6 0 1 —3 
1 -ὖ 1 —2 

0 1 -4 15 1 -1 2 

1 τ 10 

1 —2 6 — 20 1 —1 2 
αὐ 3 ue = 2 

0 1 —4 15 0 1 —3 


0 5 0 1 

0 0 1 —6 0 0 1 ete. 
0 0 0 

0 0 0 letc. 0 


Note how in each case the fanwise spread of the error 
locates its position, and that a single error grows on propa- 
gation as the binomial coefficients of increasing order. Also 
that the sum of the errors in each column crossing an error 
fan is zero. Hence, if we can establish the general trend of 
the affected differences in a column from those above and 
below an error fan, we can infer what these differences ought 
to be, i.e. what correction is needed. In the simplest case 
of approximately constant differences the errors are ob- 
tained by a comparison of mean and actual values of the 
differences. 

Thus we detect errors in a table of values from the fact 
that a single error « produces maximum errors 


€ Ze 3e θε 10¢ 20 ε 


inthe Ist 2nd 3rd 4th 5th 6th  dif- 
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ferences, respectively. This is a powerful check and should 
always be applied to any table produced by calculation. 


Example 2. Check the following table for log,~. 


ς᾽ loge differences recalculated differences 

1-0 -00000 
9531 9531 

11 °09531 — 830 — 880 
8701 133 = 8701 133 

1:2 -18232 — 697 — 697 — 29 
8004 131 8004 104 

13 26236 — 566 — 593 — 23 
7438 0 7411 81 

1.4 °33674 — 566 —512 —14 
6872 148 6899 67 

15 °40546 — 418 — 445 -18 

“ee 6454 54 

1-6 -47000 — 391 50 
6063 104 6063 | 

1:7 °*53063 — 287 — 287 — 187 
5776 —83 5776 — 83 

1:8 -b8779 — 370 — 370 177 
5406 94 5108 _ 94 

1:9 "6418 — 276 — 276 | 
5130 5130 

2:0 -69315 


The third differences indicate mistakes, probably more than 
one. Apparently 4%f,., and A%f;.g are high, 4°f;.2 and A%/;.4 low, 


suggesting too large an entry for log, 1-4; the necessary reduc- 


tion lies between 20 and 30 (as may be verified by redifferenc- 
ing), and this suggests very strongly that the digits 74 have been 
interchanged from a correct value log, 1-4=-33647. The dif- 
ferences have been recalculated after this correction, and the 
surviving error fan is re centred on 4.1; a reduction 
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of about 60 in 4fj.7 produces reasonably uniform fourth dif- 
ferences, and on checking we find that 4f;.7 should in fact be 
5716. 


Note: (i) it should be verified that the fourth differences are 
now as uniform as can be expected; (ii) a convenient check can 
be made by comparing the sum of entries in any column of 
differences with the difference between the first and last entries 
of the preceding column—this check discloses the existence but 
not the position of the mistake; (iii) differencing will uncover 
mistakes but not systematic errors. 


Problem 1. Tabulate x3— 2x°+ 3x—1 to 4D for x=0(-01)0-1. 
We could write a programme to calculate and print out the 


x f(x) "4 p= p® check 
—-02 —1-060808 — 1-0608 
30607 306 
—-01 —1:030201 — 406 — 1-0302 aif 
3020 | 302 
00 --1 — 400 —1 —4 
: θ᾽ 998 
‘01 —0-970199 — 394 ry — 0-9702 —4 
29407 | 6 294 
02 — -940792 1 —3ss“ | — -9408 «ἢ 
729019 6, 900 
03 -- -911773 — 382 — -9118 —3 
6 287 
04 — -883136 — 376 — -8831 —5 
98261 6 282 
Οὔ — +854875 — 370 — +8549 —3 
27891 6 279 
06 — 826984 — 364 — -8270 —4 
97527 6 275 
07 — +799457 — 358 — +7995 —3 
27169 6 272 
08 — -772288 — 352 — 7723 —4 
96817 6 268 
09 — 745471 — 346 — +7455 —3 
26471 6 265 
"10 — -719000 — 340 — 1190 
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value of the polynomial at each successive value of x if we were 
using a high-speed machine, but this approach is unnecessarily 
laborious if we are working by hand. For we can build up suc- 
cessive values of the polynomial by accumulation of differences 
if we start from a few calculated values and use the fact that 
7°f(x)= const. 

The section of table above the stepped line has been obtained 
by ne the terms of the polynomial, using values sym- 
metrically disposed about x=() for simpler calculation. The 
ren ‘pg is then γ᾿ ὁ" from the sequence of operations: 

7.08) = -000006, p2f(-03) = p?f(-02) + °/(-03), pf(-03) 
= pf(:02) + 7?f(-03), f(-03) = f(-02) + 7-03). This (A (8) 
dicated by arrows. 


Note: (i) that even though the third difference is considerably 
smaller than the rounding error in four decimals it cannot be 
neglected because this would constitute a systematic error which 
would accumulate ; (ii) that there is an easy check by substitu- 
tion at x=-1; (ili) that the rounded values (which make up the 
final solution) have been checked by differencing; (iv) that 
backward differences are used as these depend on values 
already calculated. 


2.6. Divided differences 

Suppose that we have to deal with a function f(x) tabu- 
lated at values a; of the argument x, where the tabular 
intervals are unequal. This calls for a more general type of 
finite difference taking account both of the change in tabu- 
lar value fti1—f; and of the length of interval x41—%;; 
the obvious choice is the ratio (fii1—f;)/(xt41—2%) which 
gives a measure of the change in tabular value per unit in- 
terval. We define the first divided difference as 


[e441] = Len) — fe) Ls), 
Mil τ Ἀὲ 
Similarly the first divided difference for any pair of tabular 
values f; and f, (not necessarily consecutive or positively 
ordered) is , 
28 
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[x,x)] ΒΟ - τὸ = [x,;,]. 
vy x i 

The second divided difference |xixt+1%1+2) of the tabular 

values f(x;), f(+141), {(x«+2) is defined as 


[χ1Δ...5] — [αι] 
gcse es ee 


χιχιιχμ5] = 
[ eexi41%1+2] ia ab 


the third divided difference as 


[¢41%140%149] — [ere %142] 
ἕω κα  ἀξνον. td ibaa Bata 


[exes X40%143] = 
X49 — X4 


and'so on for higher orders. 

The divided difference of a set of adjacent tabular values 
occupies a position in the difference table at the apex of 
the triangle of diagonals based on the tabular values, and 
its value and position are entirely independent of the order- 
ing of terms in the divided difference. This may be seen if 
the divided differences are written symmetrically in terms 
of the tabular values; for example, 


[xoxix2] 
_ [5119] — [xox] _ εἰς Ὁ —f(m)_ f(x) =~) 
Δ τα X0 X2— Kol Δ -αλὶ ΔῚ — Xo 
(0) f(x) ΠΟ 


"ζυτκυζαυ-- πο ζαι-- κο)ζαι - χρῇ. ζαῃ-- χυ)ίαε- κι) 


_ Similar results hold for differences of all orders (proof by 


induction), It can be seen from the symmetric form of these 
expressions that the ordering of the arguments is not 
Significant, and that the value of the divided difference 
depends only upon which arguments are present; thus, for 
example, [ΧΌΛΒΑΊ ΧΩ] = [xox1vo%s]. 

The nth divided differences of a polynomial of degree n 
are constant and higher differences zero; for if f(x)= 2%, 
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x" — xf 
x — Xo 


[xo] =[xox] = = x8-l 4 xM-2xy 4... + atl, 


so that the first divided difference is a polynomial of degree - 


n—1, and the result follows by repeated differencing. 


Example 3. Draw up a table of divided differences for 
f(x)=x* tabulated at x=(, 5, 9, 12, 14, 15. 


x f(x) divided differences 
0 0 
25 
5 125 14 
151 1 
9 729 26 0 
333 1 
12 1728 35 0 
508 1 
14 2744 41 
631 


18 3375 


The divided differences for a table of constant tabular 
intervai\h are simply related to the finite differences intro- 
duced in §2.3; in particular, the nth divided difference of 
n+ 1 consecutive tabular values is equal to 4*f/(n!h") (and 
it should be noted that 4*f;=7"fjin = 8"f;in,2, as these are 
alternative names for a single difference). ‘The proof is by 
induction. Suppose that [*j+1%142 ... ἀμ] = 4fi/(n! h”) 
for all 2; then 


[ 40441... Xtena] = [410442 ... Xtin+| — ΕΣ ... X40] 
(κε --- δὶ) 
= {Δ μ{π ! h*) - Anfi/(n Ι hn)} |(n 4 1h 
= An+lfi/{(n + 1)! 51}, 


The result is clearly true for first-order divided differences, 
and hence it follows that it is true generally. 


Problem 2. Evaluate (*— xo)(%— 1) ... (x—n)[*xox1%e ... Xn] 
at x= x,,1=0, 1, ... %. 
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When two of the arguments in a divided difference are equal 
the formal definition reduces to 0/0 and further consideration 
is necessary. However, 

(x— x4) ΧΧΟ ... Xp... Xn] 
= (%— 44) [4X0 ... ἀΩΔΩΔ vee XX] | 
[Ὁ ...ὄ Xt—-1% 441 ... Kui] — [NX ... Hii ... Xn] 
= (x— 1} χ-- —— 
= [x9 ...ὄ Δ Δ Χ ... Xn] — [XO ... Δ ΔΊΔΩ .νν Kn), 
and hence as x—>x;, (xn—x;)[wxo ... 4 2. Xn]—>0. 


EXERCISES ON CHAPTER TWO 


Little attempt has been made to provide straightforward exercises ον, 
this chapter as such exercises can be constructed with ease by anyone 
possessing a set of tables. Readers are urged to work such exercises so 
that they may become familiar with the various processes and the ideas 
underlying them. 

1. Show that: (i) 474=pd4f=*A, 

Gi) ΠΥ 
(Note that in each case these are merely alternative ways of denoting 
the number occupying a particular position in the difference table.) 


that: (i) A( fig) =fider+sinn ΔΙ ΔΕ δ) Ἐπ) ΔΗ, 


2 Show : i+ 
Gil) AG ie) = Afi iden) (gigi), (iv) Ad f= —Afi 8). 


; ee (—1)" n! A” 
: n in ὶ 
3. Show that: (i) 4"(1/x) ΕΝ -ὦπτεῖν 
(ii) A%e4% —(eGh —1) Max, 
where ἢ is the length of the working tabular interval. 
4. Find the polynomial of lowest degree which fits the data: 


x 0 1 2 3 4 5 
y 75 150 204 241 265 280. 
: 5. Find which tabular value is wrong in the following table and correct 
it, 
x f(x) x f(x) x x) 
1-0 "340735 1-5 ye 2:0 1: ον 
1:1 50 0 420 1°6 1293533 21 1883902 
1-2 "662671 1.7 1-435910 2:2 1-956910 
1.3 825280 1.8 1-567997 2:3 2°011438 
1-4 98609 1-9 1°687869 2*4 2046371 
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CHAPTER THREE 


Interpolation 


We assume here that tables of values provide good repre- 
sentations of well-behaved functions. Hence when we seek 
to interpolate, or estimate non-tabular values of a function, 
we take for granted that the function will behave smoothly 
between tabular points and that reasonable approximations 
can be obtained as polynomials of moderate degree. A 
polynomial of degree k will have constant kth differences; 
hence we can judge the suitability of a table for polynomial 
representation by the behaviour of its higher differences, 
and when the kth order differences are nearly constant the 
_ behaves in this range very much as a polynomial of 

eek. 


3.1. Linear interpolation 
When the first differences of a table are nearly constant 


the function which it represents must behave very nearly | 


linearly and can be approximated in each interval by a 
straight line joining adjacent tabular points; this is the basis 
of the familiar process of Linear interpolation used almost 
universally in tables of logarithms and trigonometrical 
functions. Thus in the interval x9 <x <x, 


A(x) = feo) + ——> {f(o1) —f(ao)}- 


Δ 
When we are using tables of constant interval h it will be 


convenient to write x—xo = 6(x; — xo) = 60h, where 0<0<1, 


and then in the interval (2, χε) 
(xi + Gh) = fit Ofirr—fr) = (1—- Ofc + Oca; 
or, using the various (entirely equivalent) expressions for 
the finite difference, 
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{χει 6h) = f+ OAf; = fit 67 firs ΒΞ Sit 05f;4%. 
Linear interpolation provides a simple process but one that 
can only be successful when the first differences vary very 
slowly, and this implies close tabulation to an extent that 
will be practicable only for very commonly used functions 
(e.g. ‘log’ tables). In most tables we need interpolating poly- 
nomials of higher degree and we can no longer base an 
estimate on just /wo tabular values. 


Example 4. The following extract is from Knott’s 4D log 
tables. 

ΜῈ ΝΣ. eae acta ibe A ecg τόπο κα 

10} 0000 0043 0086 0128 0170 0212 0253 0294 
11} 0414 0453 0492 0531 0569 0607 0645 0682 
12] 0792 0828 0864 0899 0934 
13} 1139 re 

proportionai parts 
8 ἢ 202 25. Fe 9 

0334 0374/4 8 12 17 21 25 29 33 37 
0719 0755] 4 8 11 15 19 23 26 30 34 


Differencing the values of the first row: 


a 6) 4 43 
100 0000 
43 
101 0043 0 
é 43 
102 0086 —1 
42 
103 0128 0 
42 
104 0170 0 
42 
105 0212 -1 
41 
106 0253 0 
41 
107 0294 -1 
40 
108 0334 etc, 


33 


INTERPOLATION 


The differences are locally almost constant and the table is 
suitable for linear interpolation: thus 


f(1033)= 0128 +- 3 x 42/10= 0141. 


The subsidiary table of proportional parts shows tenths of the 
mean difference between tabular values for a row of entries; it 
should be noted that the values printed are a compromise, and 
that the proportional parts between 100 and 102 should be 
49 13 17 22 26 30 34 39 with those between 107 and 110 as 
4 8 12 16 20 24 28 32 36, and a single table can serve the 
range (100, 110) only at the cost of some loss in accuracy (very 
nearly of the fourth decimal digit!). Good tables do not use 
averaged proportional parts because of the consequent loss in 
accuracy. 

In contrast, if we difference the first column we obtain the 


table: 

x I(x) A 43 48 
10 0000 

414 
11 0414 — 36 

378 5 
12 0792 —31 

347 6 
13 1139 —25 

322 
14 1461 


In this case we should require a polynomial of degree about 
three for interpolation, judging from the higher differences; 
hence this tabular interval is inconveniently large for a working 
4D log table. 


3.2. Interpolation using divided differences: New- 
ton’s formula 


Suppose that we have a table with constant tabular inter- 
val and we want to interpolate a value. The very act of 
putting in the ‘new tabular point’ at which we wish to 
determine f destroys the original uniformity of the tabular 
interval by producing two = intervals about the new 
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3.2 NEWTON’S FORMULA 


point, of length 0h and (1—6)h, so that we must for the 
present work in terms of divided differences. 

The problem remains essentially the same if the tabular 
intervals are generally unequal as in §2.6, and we shall 
assume this for the time being. We now seek to represent 
the value f(x) at a non-tabular value x in terms of finite 
differences for the tabular points, and this can be done in 
the following way. The original definitions for divided dif- 
ferences of successively higher order can be rewritten as 


7) = f (0) + (% — x0)[xxo], 
ΧΔῸ] = [xox1] + (x —21)[xxox1], 
[~xow1] = [xoxixe] +(% —x2)[axoxxe], 


[xxo%1 ... %n—1] = [0x1 ... a] +(*— an) [xx0%1 ... Xn). 


By successive substitution, for [xo] in the first equation 
from the second, for [xxox;] in the second from the third, 
and so on, we obtain 


f(s) = flo) + (% — xo) [xox] + (3 x0)(x— m1) [vom] + ον. 
+ (% — X0)(*— #1) ... (% τοι αὐ [Korie ... Xn] + Ra(x) 
κατ P n(x) + Rn(x), 


where 
Ra(x) = (% — x0)(* -- αι)(α — x2) ... (%— xn) [exon ... Xn). 


It should be noted that xo, x1, ... , #, can be any selected 
ordering of the n given tabular points, and that x can lie 
anywhere in the general range of the tabular points (and 
in particular need not be greater than x). 

This expression for f(x) would immediately provide an 
interpolation formula of the kind we are seeking (f(x) given 
in terms of tabular values and their divided differences) if 
only we could assert that the remainder term Rna(x) was 
small relative to f(x). When the tabulated function is a 
kth degree polynomial, divided differences beyond the kth 
are all zero and R,(x)=0 provided that n>k: hence we can 
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interpolate polynomials exactly (or, to be more precise, to 
within rounding error). 


Example 5. Find (9-3)? from the table of example 3. 

We can identify xo, x1, ... with the given tabular points in 
any order we please, and it will be sensible to select them in 
increasing order of distance from the non-tabular value 9-3 
so that the factors (x—x;) in P»(x) may be kept as small as 
possible. Hence take xp=9, x, =12, xo=5, xg= 14, *4= 15, and 
x%5=0; and then 


f(9'3)= 729+ 0-3 x 333+-0°3 x (— 2°7) x 26+0-3 x ( — 2-7) 
x 43x 1+0:3 x (—2°7)x 4:3x (—4-7)x 0+0 
= 804-357 
This result is, of course, exact (and also trivial, though it serves 
very well as a demonstration). Normally we should work to a 


specified number of decimals and there would be rounding 
error in the table and in the final result. 


The functions we have to interpolate will not usually be 
polynomials, but whenever they are well suited to poly- 
nomial approximation the higher divided differences will 
be small (and indeed this is how we judge suitability). In 
such cases we use the interpolation expression without 
remainder, 


1) =f(%0) + (ἃ -- x0)[xo%1] + (x — x0)(% — 21) [xox] + 
vee +(x — χρ)ία — 21) ... (% το Xn-1)[x0%2 ... Xp]. 


This is Newton's interpolation formula. The choice of order 
(2) for the approximation is a matter for the judgment of 
the computer and depends on the behaviour of the higher 
differences in the table. 

The truncated Newton interpolation formula 


f(x) = Ρ,( ἡ 
is a polynomial approximation of degree 2 in which the 
coefficients are given in terms of the divided differences of 
the table. The remainder R»(x) may be regarded as an 
error term; and provided that f(«) has n+1 derivatives in 
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a range including Xo, %41, ... , Xn, it can also be expressed as 


Κι(λ) Ξε (ἃ -- χο)ία τ- 21) «.. (5 τ Hn) fO(E)/(n+1)I , 


for at least one value € in the range (see problem 3 below). 
It may be seen from this form for the remainder or error 
term that for a given function we shall have a small in- 
terpolation error if the tabular intervals are all short and 
if we refrain from introducing tabular values distant from 


the level x. 


A considerable practical advantage of Newton’s formula 
is that we can at any stage add or utilize an additional value 
anywhere in the range without disturbing the terms already 
calculated. In general f(x) #Pn(x) except at the tabular 
points, but the error R,,(x) lies within assignable limits; and 
for small R, the tabular points xp, ... , x, should be chosen 


_ as nearly as possible symmetric about x. P,(x) is the 


smoothest function that agrees with f(x) at the tabular 
points, as d™+1P,(x)/dx"+1=0 for all x. 

It may be worth emphasizing that the symbol = is used 
in the expression f(x) => P,(x) to show that the right-hand 
side is a polynomial representation and that some question 
of judgment is involved in the selection of P,(x). We are 


asserting that the function is sufficiently well defined by the 
table of values and is well-represented by a polynomial of 


appropriate degree. Many authors use an ordinary equality 
sign, but we wish to draw attention to these steps of approxi- 
mation, 


Problem 3. Show that when f has n derivatives there is a 
relation of the form [xox%e ... XnJ=f( £)/n1. 

Since f(x)= P,(*)+ Rp(x) it follows that Ry has n derivatives 
too. Hence, using Rolle’s theorem, as R»(x) vanishes at the 
n+ 1 points xo, Δι, ... Xn it follows that R;(x) has at least n zeros 
in the range, R(x) at least n—1 zeros, ... and ultimately that 
RGP(*) has at least one zero, at some intermediate point ¢ of 
the range. But 

FON (x)= PLO(x)-+ R(x) =a! [on ... an] + RIM); 
and hence [xoxx2... X,]= ~ ἐν Ι. 


INTERPOLATION 


3.3. Lagrange interpolation formula 
For some applications it is more convenient to have a 
matching interpolation polynomial with coefficients which 
are given directly in terms of the tabulated values of the 
function, as in the case of the Lagrange interpolation formula, 
beige ne RO 
— (*¥—*1)(¥— He) ... (στ Hn) + 
(χο — χι)(ΧῸ — 2) ... jaar 
(*—wo)(*¥— 22)... (στο χω) 
(er —a0)(o a2) «.. (mary) * 
(x— xo)(x te δι) eee (x— Xn—-1) 
fe ee ae te Oat | Hes 
(%n—%0)(%n— Δ) ...ὄ (a — %n-ay Ὁ) 


The polynomials LZ, and Py are identically equal. For 
L,, is a polynomial of degree m which takes the values f(x;) at 
the tabular points x; (¢=0, 1, ..., 2). Moreover, Ln(x)— P(x) 
is a polynomial of degree n with the n+ 1 zeros χε, and hence | 
is identically zero;* thus Ly, is merely another form of the 


Newton polynomial Pp. 


In practice the Lagrange polynomial Ln(x) is less con- 
venient for interpolation than P,(x). The value f(x) must 
in general be more closely related to neighbouring tabular 
values f; than to distant ones; this special dependence on 


neighbouring values is recognized implicitly in the form 
of P,(x) and it is relatively easy to adjust the degree of the 


approximation to suit local conditions, while in La(x) all 


arguments appear symmetrically so that all terms must be 
calculated and adjustments are less easy. 


3.4. Interpolation in a table with equal intervals 


Most interpolation will be in tables with constant tabular 


interval and we shall now develop from Newton’s formula 


* Recall that: just one polynomial of degree less than or equal to n can 


be fitted to +1 distinct points; a polynomial of degree n which vanishes | 


at n+-1 distinct points is identically zero, 


3-4 INTERPOLATION WITH EQUAL INTERVALS 


other forms appropriate for use with constant interval 
tables. If we want the approximate value f(x) it is natural 
to choose xo as the algebraically smaller adjacent tabular 
point to x, so that x=29+ 0h, where 0<@<1 and A is the 
tabular interval; beyond this the choice of arguments need 


- not be the natural order, and indeed we can obtain a variety 
_ of interpolation formulae by taking different orderings of 
| the tabular points. 


xo Xo fo 
| Af 
Δὴ xo+h fi Afy 
4f ΓΙ 
Χῷ xot+2h = fe Af; δίῃ 
Δ | Af, 
x3 xot+3h ΚΒ Af A*fe 
3 
δ χρτ ἀδ ἢ 


" ιν - = - = * 


The simplest case is that of interpolation near the head 
of a table, in which case we identify the points xo, 44, ... of 
Newton’s formula with the tabular points in increasing 
sequence (i.e. at progressively greater distance from x). We 
must use forward differences as these provide the only 
higher-order differences available near the head of a table, 
and hence x,=%+hkh, x—x,=(0—k)h. Bearing in mind 
that [xox ... %n]=A"fo/(m! h"), we have by direct sub- 
stitution in Newton’s formula 

.θ θ(θ -- 1) ., θ(θ --- 1)({0 -- 2) 
Κῶ ΞΡ τι 4fo+ ---οτ-- A*fo+ γεν χ re A8fo + 
4 (9-1) .-» (0 --- 
᾽ν (n+1)! 


This is the Newton-Gregory forward difference interpolation 
formula for use near the start of a table, mainly in the 
interval (ao, Δ1). 
In a similar way we can construct a Newton-Gregory 
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extending a table beyond its 
*=Xn_+ 6h, we obtain 
F(x) = fnt+ ὅν, + At θ) »ϑ, + woe + 


" θ(1-Ὲ θ)... (κα -- 1 - θ) 
n! } 


f(xn). Thus, 


V"fn. 


Most interpolation will be in the interior of tables, where. 
we use central differences so as to give equal weight to the 
tabular values on either side of the working level. In this 
case we can introduce tabular values roughly symmetrically 
about the interpolation level if we identify the xo, x1, ... of 
Newton’s formula with the tabular points in order of 
increasing distance from x starting with the point alge- 
braically just below and alternating above and below, as 


“o-3h fs 
x4 wo-2h fis sy 
XQ wo-h fy > oe 
df-1/2 δ... 
x9 Xo fo δὲ d4fp 
δῇ, δϑῇ 2 
δὴ woth ἢ δ Sf, - δέῃ 
x3 to+2h fo = δὲ fan 
df5,2 
x5 *o+3h fy 
xo+ 4h ta 


illustrated. In this case δ τα Xo + Oh, *—xo4=(0+R)h, 


* Extrapolation is a potentially more dangerous pursuit than interpola- 
tion and we can relatively easily extrapolate functions into regions wind | 


they do not exist! 
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x—Xer-1=(0—k)h, and — [xox ... xen] = 52%fo/(2n)! h2, 
[ΧΟΛῚ ... Xon41] = 62" +1 fiy/(2n +1)! h2"+1 on account of the 
selection order of the tabular points; and on substituting 
in Newton’s formula 


f(x) = fot (*—xo)[xo%1] + (α —x0)(* — 21) [xoxx2] + ... 


we obtain the ‘basic Newton formula’ for constant interval 
interpolation 


T(x) = fo O5f44 CD δέ: OK 1.611) δϑβμ + 
OEE NOH?) sty 


Note that if we sometimes write our formulae as unter- 
minated series this does not mean that they are to be 
regarded in any sense as infinite series, but only that the 
truncation point is left to the discretion of the computer in 
each case; if polynomial representation is suitable, the error 
will at first decrease as we take extra terms, but beyond a 
certain stage it will generally increase again. The basic 
Newton formula can be rearranged in more useful forms. 
(1) For interpolation over the inner part of the interval 
xox, we can derive a form symmetrical about the mid- 
point x9+4h by rewriting even differences as, 6.5. 


B2fo= 4(B4fo+ BY) + 4(3%fo— B%f1) = H(5Yo + 84/1) -- φδϑι. 
Thus 
fx) = U(fo+-fi) +4 fo—fi)} + O8fin+ 
+ ED (4820+ 95) +484 BY) +... 
A(fo+fi)+(0-H8fu+ SE) (δε 9%) + 
ε{0- 8 - Ὁ sy, 
, (6+1)0(0— 
* 241 
D 41 
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This is Bessel’s interpolation formula. Bessel’s interpolation 
formula is sometimes written 


His) = fin (0-4) 8fn τι Ὅμδο, + EDOM Τρ, 
+ ose 
where the averaging symbol » denotes the mean of two 


adjacent elements of a vertical column; thus 2ufa=fo+fi, 
etc. It is also written 


70) = fo+ O8fta + B'(S%fo+ 8/1) + B"H fin + BY(S4fo+ S4f) 
+ ass 


erp 


where B"=46(0—1)/2!, B” = 6(0—4)(@—1)/31, ... are the 
Bessel coefficients. 
(2) Some sets of tables include auxiliary sub-tables of even- 


order central differences (as these lie on tabular lines); to _ 


obtain an appropriate interpolation formula we need only 
replace the odd-order differences in the basic Newton 
formula by the one-lower-order even-difference combin- 
ations, 


Hea) = for (fr— fo) + Hoy 5 65 NON 504, — sap 
T see 


Taking $= 1-0, 
Hee) = $4 for =F * sxf04§ FPP oy, 
θ: 


εὐ fe δι CoE ap sd 


This is Everett’s interpolation formula. 
(3) We can obtain a form symmetrical about tabular 
points as follows: 
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f(x) =fo-+ O(8fm—48%o) + 5 0°B4fo-+ 4 0( 82 1)(5%fu— 3 δι 0) 
+ 2 62(02— 1)54fy+ ... 
= fot 40(5f-m + Bf) +3 O3%fo+ 
* xa 1)(53f_m% + S8fth) + ..Ψ. 
=fot+ OpBfo + 7 0°B%fo-+ 700 1) 15%fy + 
+ 20402 1)8Yfot su 


and this is Stirling’s interpolation formula; we shall use this 
form later in deriving differentiation and integration pro- 
cesses. 

A useful formula for half-way interpolation follows from 
the Bessel or Everett forms by taking 0=4: 


Fin 2 3 (for fi) g(O%fo+ BY) + ;59( 84+ BYfi)— 
— Faq H+ δ} + OC, 


This is simple to apply and useful as a first step in breaking 
down a table with large interval. Note that this form of 
expression indicates approximation and also shows (in 
terms of the interval) the order of magnitude of the error 
term. 


Example 6. Use Bessel’s interpolation formula to find 
f(10/3) from the table for f(x)=1/x given on page 19. 

Here 6=(10/3—3-3)/(3-4—3-3)= 1/3, so that the coefficients of 

the Bessel formula are 6(0—1)/2.2!= —1/18, 0(8—4)0—1)/31 
= 1/162, (8+ 1)0(@—1)(8—2 ae l= 5/486, ... Thus, 
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10... 1 1 
ΑΚ 3 ) = -303030+ 3( —°008912)— is ‘000558 + -000508) + 


1 | 
Ἔτροί —*000050)+ ... 


= *303030— -0029707—-000059e— ‘0000003+- ... 
τ -2, 


i.e. {(34) = -300000 οῃ rounding off, 


Note: (i) that Bessel’s formula up t di is j 
; ' Ρ to first differences 
linear interpolation, and is most easily calculated as Sot ἜΝ 
(ii) the guard figure kept as a subscript to minimize rounding 
error (which might otherwise accumulate undesirably) ; (iii) that 


re a from 8%f is just worth retaining, that from 54f 


3.5. The use of 
throwback 


For an isolated interpolation we shall 
with which we are most familiar, 
polations 
formula, 

Bessel’s formula is 


Bessel and Everett interpolation: 


use the formula 
but for a run of inter- 


working 
of dif- 


finding the greatest magnitude of each Bessel : 
0<9<1; thus max|6(@~ 1)/21] =-125, essel coefficient for 


ence we | 
neglect S<4, 33f<60, 54f<20 units in the last dedienal 
with consequent error less than a 4 unit in the last decimal 
error in any particular 

can be found from the error term, but shank: seegti valtios 


it will be advantageous to select the ‘best? — 


| 
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BESSEL AND EVERETT INTERPOLATION 


If only first differences are given the table will take linear 
interpolation ; but when higher differences are needed it is 
common to print only δὲ and 54f because of the difficulty 


involved in printing odd differences on half-lines, and in 


this case Everett’s formula will obviously be more con- 


_yenient. A further advantage is that second difference 


terms in Everett’s formula take complete account of 8%, 
and similarly Everett up to fourth differences actually in- 
corporates the effects of the first five differences; hence it 
will be superior when third differences cannot be neglected 
in Bessel’s formula (5¢f>60) but fourth differences are 


negligible, or when fourth and higher differences are sig- 


nificant (as when many decimal figures are needed, or 
when the tabular interval is large). It can easily be shown 
that δὲς 4 and 8%f<20 can be neglected in applying 
Everett’s formula with a consequent error of less than half 
a unit. 

In an effective interpolation formula we hope that there 
will be few terms to calculate and yet we must avoid exces- 
sive truncation error, and we can make gains in both these 
directions by taking partial account of one of the higher 
differences in the procedure known as throwback. For 
illustration, in Bessel’s formula the whole contribution 
from the second and fourth differences is 


OF (Yo 8%/i)+ (δι: δι) 


= TAT [δ τ BA) -- γχα + θχ2-- OV oY + δι). 


(0+ 1)0(9-- 1)(9-- 2) 
2.4! 


By a fortunate numerical chance the coefficient 
(1+ 0)(2—6)/12 varies little over the range 0<@<1; in- 
deed it lies between -1800 and -1875 for -2<0<-8, and 
near the ends of the 6-range the whole contribution is 
small because of the outer factor 6(@—1)/4. Hence we get 
& very good approximation ay Bessel formula carried to 
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e(1+9)(2-9) 


fourth differences if we work to third differences in the — 


formula and use the modified second differences 
ont {:- Bf; δὰ Corp, 


where (ὦ is a constant, and may be taken as 0-184; this is 
known as throwing back the fourth difference on the 
second. Some accuracy is lost (though this throwback error 
is <4 unit for 64f=1,000 units!), but we make a sub- 
stantial saving in labour, and this is precisely the sort of net 
balance we always seck in designing numerical processes. 

In the same way we can throw back fourth differences on 
the second in Everett’s formula, and again we use 3,2/; 
= 82/;— 0-184 54f; in place of 82f; and the modified Everett 
formula is 
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f(x) = (1—A)fot+ θ6ῃ -- 1θ(θ —1)(6 — 2) 82, fo+ 
+4(0+ 1)0(0—1)82, fi. 


Summarizing: (i) near the start or finish of a table use a 
Gregory formula; (ii) when 6°f<60 units in the last digit 
use Bessel up to second differences; (iii) when 6°f>60 but 
54f<1,000 units (or 400 units if 54 and 544 are of opposite 
sign) use Everett with throwback from fourth to second 
differences; (iv) when 54f>1,000 units use throwback with 
higher differences. 


- 7. Find (123) from the table for f=1i/x given 
ow. 

In the following table h=-05 and x= 2/3=-65+h/3, so that 
6=1/3. Third differences cannot be neglected, but as fourth 
differences are less than a thousand we should normally use 
Everett with throwback, and we need then take no further 
account of δέ. As an illustration this particular example will be 
worked using : (i) Bessel to 54f; (ii) Bessel to δ3 with throwback ; 
(iii) Everett to 82, f. 


x f(x) ὃ δὲ ὃ9 δέ |.-- (δ δῆ, ομδϑ, 


60. 2-00000 
— 18182 

δῦ 181818 3031 
— 151651 —701 

60 1-66667 2330 203 | —37 2293 
— 12821 — 498 

65 1-53846 1832 131 | —24 1808 
— 10989 —367 

‘20 1:42857 1465 93 | —17 1448 
— 9524 —274 

‘75 1:33833 1191 63 | —12 1179 
— 8333 — 211 

80 1-25000 980 
— 7353 

"85 117647 
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Hence: (i) 
“O-1)/4=  —1/18 
6—4)(0—1)/6= 1/162 
(0+ 1)6(0—1)(0—2)/48= 5/486 


O8fi,= —+036630 

— (δὲ -Ὁ &f,)/18= —-001832 
53f1,/162= —-000023 
5(S4/o+ 54f1)/486=  -000028 


1-49999s 
(i c) 
43-- 6 —1/81 
hae en 7 
fo= 1-53846 fo= 158846 
SfuJ3= —-036630 Sfv,/3= —-036630 


—(8,2fo+5,3/1)/18= —-001809 —5,2fo/81= —-001116 
Care alae — 000023 — 45,2/:/81= —-00071s 


1:49999s 1-4999% 


In each case the rounded solution is 150000. 


Note: (i) that the subscript m denotes a difference modified 
by throwback; (ii) that the effect of neglecting δῆ is about 2 in 
the last place, and some account must be taken of df if the 
interpolate is to be as accurate as the tabular entries; (ili) that 
successive terms do not necessarily decrease monotonically in 
magnitude; (iv) that the probable round-off fluctuation of +6 
in 54f has negligible effect. 


3.6. Some uses of interpolation 

(i) Subtabulation. The choice of tabular interval in the 
construction of a table is influenced by the uses anticipated 
for the table, but also by the economics of book production; 
hence we shall sometimes find that the printed tabular 
values are too widely spaced in regions of special interest. 
The process of ‘filling in’ between tabular values (known 
in this context as pivotal πιο by systematic interpolation 
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is Called subtabulation. For example, if the initial tabular 
interval h is to be broken down to 4/10, it may be worth 
while applying the halfway interpolation formula followed 
by Everett using throwback as necessary. Subtabulation 
should be carried at least to an interval such that accurate 
interpolation is obtained using no higher than second dif- 
ferences; once a table can be linearly interpolated there is 
no gain from further subtabulation. 


Note: (i) that if a table of a complicated function is to be 
built up by calculation from a formula and a fast computing 
machine 1s not available, the easiest procedure may be to cal- 
culate at large interval and subtabulate; (ii) that after halfway 
interpolation the differences in a table will be quite different 
and higher differences will be considerably reduced (as 2-"), 
and that the usefulness of subtabulation arises from the fact 
that 6"f=((A"); (iii) that a table of irregularly spaced values 
can first be reduced to regular interval using Newton’s divided 
difference formula, and then subtabulated in the normal way. 


(11) Inverse interpolation. In the process of interpolation 
we seek the value of f(x) for a specified non-tabular value x; 
the inverse process in which we seek the value of x cor- 
responding to a specified non-tabular value f(x) is also im- 
portant (e.g. in extracting roots of the equation f(x)=0), 
and this is called inverse interpolation. It can be carried out 
by successively closer subtabulation in progressively 
smaller neighbourhoods of the specified value f(x) until 
linear interpolation is possible (note that this approach 
allows full scope to the computer for an inspired guess at 
the neighbourhood and interval for starting). 

Inverse interpolation can also be carried out by methods 
of successive approximation; this approach can be illus- 
trated using Everett’s formula. Everett’s formula to second 
differences can be rearranged in the form 
6= ἰὼ —for+ Am )@- 2)8*fo — (0 + 1)82fq] + sh 


1 . 
fi-f 
Provided that the function f(x) is well defined by its table 
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of values, the terms involving second and higher differences 
will normally decrease fairly rapidly, so that we can expect 
to get a modest first approximation @, by ignoring the con- 
tributions of second and higher differences, and taking 

9, f2)—fe _ flx)—fe 

aie x. δβμ "ἢ 

To improve this approximation we must bring in the term 
in second differences, and this involves @; but as this term 
will provide only a relatively small correction to 6, we can 
get a satisfactory second approximation 62 from 


Oa = {f (2) -- [0 + ἐθιίθι — 1)[(01 — 2)8%fo— (θι + 1)δ 51} |δύβι 
where the error introduced by writing @; in place of @ in the 
second-difference term is likely to be of the same order as 
the higher-difference terms which are neglected. The 
accuracy of the approximation can be improved by the 
successive introduction of terms in higher-order differences 
(though only to the degree permitted by the quality of the 
table). Note that the number of decimals which should be 
carried at each stage can be estimated by making the round- 
ing error roughly comparable with the first term neglected 
(and it should be stressed that there is no gain from carry- 
ing too many figures in the early stages of the calculation). 


Example 8. Find the value of x for which log, x=0°3 from 
the following table. 


x log x ὃ δὲ 58 δ΄ δδ δ 

12. 182822 — 6970 — 292 — 57 
80042 1036 rons 

18 *262364 — 5934 — 217 — 16 
74108 819 59 

14 +336472 — 5115 — 158 — 23 
68993 661 36 

15 405465 — 4454 — 122 —3 
64539 539 33 

16 "4τ0004 — 3915 — 89 —21 
60624 450 12 

1:7 +530628 — 3465 mitt +10 
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Working from the table, 

ὃ = O3—fo__—*3— 262364 
fi-—fo  *336472—-262364 
Our best choice is probably to proceed with the value @,=0°5, 
if we take due account of ease of calculation at the next step and 
also the likely error in this term! Then, using Everett up to 
modified second differences (we may just as well introduce at 
this stage the small but not negligible contribution from the 
fourth difference) we have for the second approximation, 


θ5-- { f(x) — fo+ 401(61— 1)[(61— 2)832,fo— (01+ 1)82,fi]}/Sfin 
-(3- 262364— 1|- - -005894)— s(- 005086)|| + 
“074108 


= '50785s. 


= 498593. 


A recalculation of @2 using the improved approximation 

§<=-49859 in place of 0<-5 will now serve both to reassure us 

that no substantial error has arisen from use of the crude 

value ;=-5 in calculating θᾳ (i.e. that the approximation is in 

το ry eee to the second stage), and to check 62 against 
unders. 


6,= {-3—+262364-+ ἔχ -49859 x (—-50141)[—1-50141 x 
eed x (— -005894)— 1-49859 x (— -005086)]}/-074108 


There will be no gain in accuracy by proceeding further, as the 
fourth differences (of about 200) have been sufficiently taken 
into account by throwback and the contribution to the numera 
tor from the sixth differences is negligible. Hence 


# =1-3+°49859 x 0°1= 1:349859. 


_ Note: (i) that the irregular appearance of the sixth differences 
is sufficiently explained by the probable round-off fluctuation 
of +22; (ii) that by recalculating 02 we avoid carrying on an 
unnecessary error to the next stage, and in this case if we had 
used #,;=-508 we would have found #2=-49859 and might have 
carried on -49860; (iii) that processes of successive approxima- 
tion in which each improved estimate is based on the previous 
one are called iterative sii ᾿ (iv) that modified differences 
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are useful in inverse as in direct interpolation; (v) that the 
accuracy of the process is good provided that δύῃ is not small, 
but that for small 8/,, the rounding errors can be disastrously 
amplified, for example near a stationary point of f(x), and it is 
then usually best to subtabulate until ἤρου: interpolation is 
possible ; (vi) that an initial survey of the problem is absolutely 
essential so that the computer may have some idea of the sort 
of pitfalls that must be avoided ; (vii) that when higher differences 
are large it may be best to subtabulate at interval 4% as a first 
step; (vill) that on no account should the interpolation formula 
be treated as a quadratic (or cubic) to be solved for @, as this is 
a laborious and inaccurate approach. 


3.7. General notes on interpolation 


(i) Tables of Bessel and Everett coefficients for @2=0(-001)1 
can be found in Interpolation and Allied Tables (H.M.S.O.), 
in Chambers’s Six Figure Mathematical Tables, Vol. 2, and 
elsewhere. 

Critical tables of coefficients are specially useful, and a 
set of critical tables for Bessel coefficients is given in 
Chambers’s Shorter Six Figure Mathematical Tables. Bessel 
coefficients vary slowly with 0, and it is possible in quite a 
short table to write down ail the values taken by a co- 
efficient, for example B”=-001, -002, -003, ..., and to 
list beside them the critical values defining sub-intervals 
of @ which correspond to each value of Β΄, The leading 
terms in the critical table for Β΄" are: 


θ Β" 
Ὃ2 ὀ- 000 
— 02 
ἘΠ 003 
‘O14 rt 


These critical tables are easy to work from and are very 
convenient. 
(ii) If instead of a table of values we start from a function 
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which can be differentiated with reasonable ease it may be 
more effective to evaluate the function and its first few 
derivatives at pivotal points of the range, and then to sub- 
tabulate using ‘Taylor series based in turn on each successive 
pivotal point. The terms in the Taylor series are in general 
smaller than corresponding terms in any finite difference 
interpolation formula. , 

(iii) Scientific literature contains many accounts of (so- 
called) semt-empirical theories, in which dimensional and 
other incomplete arguments are used to derive solution 
formulae containing several disposable constants, and these 
constants are to be evaluated from experimental results to 
complete the solution, But amy given segment of curve can 
be matched to specified accuracy by a polynomial, which is 
fully determined by its coefficients. Put in another way, any 
curve segment (or set of experimental points) can be matched 
by a function containing a sufficient number of disposable 
constants, and the relatively simple curves which ought to 
be the product of well-conceived experiments can normally 
be matched over moderate ranges (and within experimental 
accuracy) by quite low degree polynomials. Hence we must 
treat semi-empirical theories with great reserve; for the 
link between theory and experiment has been reduced (at 
worst) to one of curve-matching and we have forfeited the 
direct relationship which alone justifies our construction of 
theories. 


EXERCISES ON CHAPTER THREE 


Again exercises may easily be constructed from any set of tables. 
1. Subtabulate the following table to produce a table of f(x) for 
*=1-5(1/30)2-4. 


x de x F(x) x (x) 
1.5 546 1.9 64185 23 3291 
1:6 «47000 20 “69315 2-4 “87547 
1:7 53063 2"1 .74194 
1.8 58779 2:2 -_ 
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2. State which of the two following tables is suitable for linear mter- 
polation and estimate the subtabulation needed in the other to make it 


suitable, 
(i) "176809 (ii) “176916 


EXERCISES it 


7. Show that the remainder term for Bessel’s interpolation formula 
can be written as 


177840 177896 6(@2—12)(62—2) ... (@—n—1 : if 
178872 178891 OPE) oe On) han f@n) (ὦ), | i 
179904 179904 | 
*180937 180937 and find the corresponding remainder terms for the Everett and Stirling | 
"181970 "181993 interpolation formulae. | 
183003 183075 
184037 184186 
3. Find f(4-2) from the table. 

x f(x) x F(x) x Fix) 

20 = *500000 40 +250000 60 166667 

2:5 “4.00000 45  +222222 65 153846 

$0 +333333 50  +200000 

35 = 285714 55 181818 


4. The following extract is taken from the table of natural tangents 
in Castle’s Five-figure logarithms and other tables. 


: i ae Nae AS Raa MN il aia 
67| 2-35585 36733 37891 39058 40235 41421 42618 43825 45043 46270 
68] 2-47509 48758 50018 51289 52571 53865 56170 56487 57815 59156 


Ὶ mean differences 

Ἧς ie 2 3 4 5: 

67 199 397 596 795 994 

168 mean differences cease to be suf- 
ficiently accurate 


Find whether linear interpolation does in fact provide 5D accuracy, and 
determine the errors which will be introduced if the mean di 1088 
are used in calculating the tangents of 67° 3’, 67° 27’ and nd 67° 57’. 


5, Evaluate tan+/x when x=-075 radian, using the table for tan +/x 
_ given in problem 4. 


6. Subtabulate the following table to produce a table of f(x) for 


x=0(-1)2°5. 
x x) x f(x) x F(x) 
0 1 1-8 3°10747 23 503722 
0-9 1-43309 1.9 3.41773 2.4 5°55695 
1-4 215090 21 4-14431 
16 2-57746 22 456791 
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Differentiation and Integration (Quadrature) 


The process of interpolation consists essentially in fitting 
a polynomial to neighbouring tabular values of the func- 
tion, and then in substituting the interpolation polynomial 
for the known or unknown function which the table of 
values is supposed to represent. ‘The error of approxima- 
tion involved in using P,(«) in place of the actual function 
f(x) can be kept within any prescribed limits over the whole 
of a (finite) x-range* if we use polynomials of sufficiently 
high degree, which in practice implies retaining a sufficient 
number of terms in an interpolation formula. The general 
character of this process is illustrated in figure 2 where the 
continuous curve, which represents an interpolation poly- 
nomial P(x), oscillates about the dashed curve representing 
the true function f(x) and intersects it at tabular values 
(marked by crosses). The lateral separation has been mag- 
nified so that the figure may be clearer. If we wish to form 
derivatives or integrals from a table of values we have no 
option but to base our numerical processes on an approxi- 
mating polynomial P;(x) (or some other approximation); 
otherwise we cannot find the derivative, for example, even 
at a tabular point. 

Although the difference between P,(x) and f(x) will be 
small over the working range of the 5 epee and 
zero at the tabular points, the gradients P,(x) and f’(x) are 
likely to differ more markedly, particularly near tabular 
points. This effect can be illustrated very simply by the 


4 Thus the approximation is uniformly valid over the range. 
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non-polynomial (and in fact quite useless) approximation 
S(x) =f(x) + ε sin(27x/h) to f(x), where ε is small and his the 
tabular interval. The greatest difference |.S(x)—f(x)| maz 
between /(x) and its approximating function is ε at each 
ialfway point. However the greatest difference between 
the gradients |S’(x)—f'(x)] ma 5 2zeih) cos (27x/h)| maz 
is 27e/h at tabular points, and between Ath derivatives 
| S(x)—f(x)| maz is (27/h)*e at tabular points for 
odd-order and halfway points for even-order derivatives. 
Thus for h=-1 the approximate first derivative is in error 
by about ὅθε, and the Ath derivative by about (60)*e |! 
Thus numerical differentiation is a process in which 
accuracy is Jost, and in ee the determination of 
κ 
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higher derivatives numerically is usually unsatisfactory 
and is to be avoided if any alternative process is available 
(as by successive differentiation of the function before 
insertion of numerical values, or by successive differentia- 
tion of a differential equation, so that higher derivatives 
may be calculated directly by inserting numerical values of 
quite low order derivatives, etc.). 


ὃ 
The numerical evaluation of an integral | f(x)dx is often 


termed quadrature to distinguish it from the integration of 
differential equations. The integral over any interval can 
be interpreted as the area under the curve in this interval, 
and it follows that integration involves an averaging of the 
error f(x) —Pn(x) in each tabular sub-interval and very fre- 
quently substantial cancelling of error contributions from 
adjacent sub-intervals (as within the working range of a 
polynomial approximation the error is alternatively + ve 
and —ve between successive tabular points), When we 
accumulate contributions to an integral over many sub- 
intervals the overall accuracy will be greatly enhanced 
because of alternations in sign of local error terms. Thus 
quadrature is a process in which accuracy is gained (in 
powerful contrast with numerical differentiation), and this 
effect is sufficiently marked in practice for us to adopt a 
basic integration process that depends on quadratic (i.e. 
parabola) matching. 


4.1. Numerical differentiation 


Formulae for numerical differentiation can be obtained 
easily by differentiating the interpolation polynomial, and 
this can be carried out conveniently using one of our inter- 
polation formulae. Derivatives are needed most commonly 
at a tabular point, and hence from §3.4 we select Stirling’s 
interpolation formula as the form symmetric about tabular 
points. When this is a 


ΤΩ =) = 7 | Bf + 65%fo+ < | αὐλῶνας 


HOP Davy. 


f"(*) = aan Sfp ok ude + Aor athe ΒΕ 


PED eal 


and so on for higher derivatives, where x=xo+ 0h. These 
forms can be used directly to obtain derivatives at non- 
tabular points by suitable choice of 6; and with @2=0 we 
obtain finite difference formulae for derivatives at tabular 
points as, 


hfo'==pdfo - ἐμδϑῇ + 35 μδϑῇῃ Ὁ- Tpit + sees 
hfe" =8%fo— T5!fo+ pith Ἄ ε: μὰ... 


It must be emphasized again that these are not infinite 
series, but that the computer must select the truncation 
point appropriate to each application. Thus if a poly- 
nomial representation of degree five is appropriate for 
some application we use the form 


ye 1 1 ) 
fo =; (udp — gH δον sgh 
which can also be written in the form:* 
Tee ener 
fo =; {Sf geist 3 p18} + OC). 
The advantage of the latter form is that it gives an im- 


5. This form implies that in the case of a function suited to polynomial 
representation the first term neglected (i.e. 46"f,/140h here) is of order 
ἐδ, and gives a good estimate of the total error for small ἢ, 
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mediate indication of the effect of a change in working 
interval on the truncation error. 

It should be noted that these differentiation formulae 
are of the general form f(")= (some combination of dif- 
ferences)/h*, The rounding errors in the (combination of 
differences) depend very little on the interval ἡ, so that the 
consequent error in the derivative is proportional to A-*, 
and will be greatly amplified even for low-order derivatives 
when ἡ is small. Hence in differentiation we should use as 
large an interval as is consistent with a proper representation 
of the function, even though this may involve the introduc- 
tion of a number of higher-order difference terms. 


Example 9. Find d(cos x)/dx at x=0-5 radians from a 6D 
table of cos x with interval 0-001. 

From our experience with interpolation we might choose to 
work with as small a tabular interval as possible, -001 in this 
case, Extracting that part of the difference table actually needed, 


δ: F(x) δῇ Sf 
‘499 -878062 -1 
- 479 
00 υ877ῦ88 -1 
— 480 
01 877103 1. 


The general formula for the first derivative is 


Ife pdfo— =udYot aap B5fo— ann 


and as third differences are negligible we shall truncate it after 
the first term, so that on substitution from the table -001 f’(-5) 
= — ‘0004795, and f’(-5) == —-480. Note that the maximum 
error in #df from rounding of the tabular values f(x) is half a 
unit in the sixth decimal, and that in consequence f'(-5) has a 
possible error of half a unit in the third decimal place! 

Repeating the calculation with tabular entries at interval -01 
selected from the table, 
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‘49 = 886995 — 88 
— 4750 

‘50 °877583 —88 
— 4838 

‘Sl = *872745 — 88, 


we again use the formula Afo=~pdfp and hence -(1f'(-5)= 
—°004794, and f’(-5)<=—-4794 with a possible round-off error 
of half a unit in the fourth decimal digit. Again with h=0-1, 


"Δ 921061 — 9203 93 

— 43478 434 —5 
ὃ "877588 — 8769 88 

— 52247 522 —4 
‘6 °825336 — 8247 84 : 


“1f'(-5)+— -047862s— -0000797— 0000002, and f’(-5)=—-47942 
with a possible error approaching one unit in the fifth decimal 
cae from the round-off error in the two terms of the formula 
ify = pdfo— 42.55/o). Finally with h=0-3, 


‘2 980067 — 87547 7822 
— 102484 9155 — 820 

‘5 *877583 — 78392 7002 — 622 
— 180876 16157 — 1442 

‘8 -696707 — 62235 5560 


3f'(-5)==—-141680—-0021093—-000037s+ ... . Still further dif- 
ferences must be taken into account (i.e. a higher degree match- 


ing is needed), and the difficulty is that at the term aap B8fy 


the contribution of the probable rounding error is about 4x 10-? 
and is increasing more rapidly than the coefficients of the series 
decrease; hence we have reached a stage at which the rate of 
decrease of higher-order differences is too slow for numerical 
differentiation (i.e. the tabular values utilized do not define the 
function sufficiently well for this purpose). 


Note; (i) that the coefficients of the differentiation series 

decrease rather slowly; (ii) that there are two distinct dis- 

advantages of numerical differentiation, the amplification of 

round-off when the tabular interval is small and the succes- 

sively poorer matching of the a P(®)(x) and f{()(x) as 
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the order & increases; (iii) that the calculation of f’(-5) from the 
above data goes much worse, (-(001)%f"(-5) == —-000001, 
(01)*f"(-5) == —-000088, etc. 


4.2. The integration of a specified integrand; Simp- 
son’s rule 


The process of quadrature, that is the numerical evalua- 
: | b 
tion of the definite integral | f(x)dx of some specified 


εἰ ὦ 
function f(x), is of considerable importance, as we are often 
unable to find effective analytical processes for evaluating 
integrals even though the functional form of the integrand 
is known. The obvious way out of this difficulty is to ap- 
proximate to the known integrand f(x) with some function 
we can integrate, and we are immediately on familiar 
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ground as it was for a very similar kind of approximation 
(to the supposed function underlying a set of tabular values) 
that we have introduced interpolation polynomials, and we 
are really rather good at integrating polynomials! 
The simplest of interpolation polynomials is the chord 
joining two points of the curve y=/(x), 
y = Pi(x) =fo+ O4fo. 


Hence we have the crude approximation 
ath a th 1 
{ f (a)dx=| Ρ i(a)dx=h| (fot+ 0Afo)dé, 
* a 0 
or Ἷ 
2+ 
[Ὁ Heya io+fd, 
af) 


which is known as the trapezoidal rule, and consists in re- 
placing the integral by the area of the trapezium ABCD, 
as shown in figure 3. A substantial interval a<x<b can 
be subdivided into a number of sub-intervals of length ἡ 
in each of which the simple trapezoidal rule is applied, 
to give 


bez +nh } 
oN γο)άχ ἘΠ ἘΠῚ (ον ἘΠ  Ἐ3ὼ. 
= 0 
The trapezoidal rule will obviously give a poor approxima- 
tion unless the sub-interval ἢ is short enough to identify 
each time an almost linear region of the curve y= f(x). 

A much better approximation can be obtained by fitting 
an arc of a parabola to the curve y=/f(x) at the three 
tabular points xo, xo +h, x9 + 2h; that is by using the inter- 
polation polynomial P2(x) in an interval of length 2h. Then 


ao 
=| “ot Poa) = hy \fo + OAfo+ ah a6, 
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%t2h russ 

[Ae ydea Gio Hit fd. 

This is Simpson’s Rule, and it provides a remarkably simple 
and effective process for numerical integration. We shall 
show that the direct application of Simpson’s rule will in 
many cases provide sufficient accuracy, provided that the 
working interval A is sufficiently small; this will usually 
involve subdivision of the whole interval (a, b) of integra- 
tion into an even number of sub-intervals, when 


[ ΠΝ fede SH for 4fs + 2fa+4fs+ fat... +4fon—1+fon). 


A variety of quadrature formulae can be obtained by 
using interpolation polynomials of successively higher 
degree with matching at a successively larger number of 
tabular points, or by suitable combination of other formulae 
in such a way as partially to cancel errors. However, these 
tend to be very largely of theoretical interest, as they in- 
volve subdivision of the range into awkward numbers of 
intervals (which are integer multiples of the degree π of the 
matching polynomial P,(«) ), and also they tend to involve 
a variety of coefficients so that they are less convenient in 
use and more liable to blunders. For example, the three- 
eighths rule is obtained using the matching polynomial 
P,(x) fitted at xo, xo+A, x0+2h, xo+3h; and Weddle’s 
rule using Pe(x) fitted at xo, ..., *o+6h (see exercise 3 
for these formulae). 

We must have an estimate of the error involved in these 
quadrature processes before we can apply them confidently. 
Suppose that f(x) is continuous in the range (xp —h, x0 +A), 
and in this range has the Taylor expansion 


fe) = fo afo! + pao" + aif + wu. + Row 
where Roz =x2*f(2h)( oe some £= é(x) between 


4.2 SIMPSON’S RULE 
xo and x. Then on integrating over the symmetrically dis- 
posed range (xo -- ὦ, xo+A), 


ath 2: 1}2 ἐφ 1 hs ὥ 1 Soe i 
J? fovie~ 2M έτη ote a O°"? + aes ay } “Ἀμάν. 


Simpson’s rule can be μὰς Spas over the same range in 
terms of the derivatives f{ by expanding fi=/f(x%o+h) and 
fa=f(xo—h) separately in Taylor series and combining 
the two terminated series, 


[θα = i f-a+4foth) 


τρίτης os ta 


3 (aay FOMCE) 
+ Yor forfeit Lefe+ . χω 8) 


| Lh? . 1.8. 
= ; = =f, — —-fal4) 


where xo—h< £:<% and x0 < &2<%0+h. By comparing 
these two forms it may be seen that the correction needed to 
Simpson’s rule to give correct quadrature is 


eit 5 Yur Ble 7f,(6) — 79"? 
5p fo" i300" ‘fo — ... + ᾿ Cyne ih 
o 


- πῆρ i a+1ff (2%) ( £1) + f (2 k)( £2)}. 


If f@*)(x) is bounded in (x9—h, xo+h) it is easy to show 
that the remainder term in this correction has the value 


—F(k—1)h#+1f2%)( £)/(2k+ 1)! for some € lying in the range 
(xo—h, xo+A) but not sie equal to any of the 
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previous values. Provided that the interval h is chosen so 
that f(x) is well-represented by a table of values at inter- 
val h, the first of these terms will dominate the correction 
and yet will itself be small. Thus we can integrate a very 
wide range of functions with Simpson’s rule provided only 
that we are prepared to use a suitably chosen working 
interval, and in most cases this will involve subdivision of 
the range of integration; the extent of this subdivision 
will be a matter for judgment based on the character of the 
integrand, We shall also define the error = —(correction). 

The accuracy of formulae based on higher degree inter- 
polation polynomials can be estimated in a similar way (see 
exercise 3). Formulae based on polynomials of even degree 
are found in general to be more effective than those of odd 
degree; thus the corrections expressed as single terms, are 
—hf,(€)/12 for the trapezoidal rule, —h5f(¢)/90 for 
Simpson’s rule, and — 3h5f{1)( £)/80 for the three-eighths 
rule (where the ranges of integration have been taken as 
h, 2h, and 3h, respectively, and in each case & is some point 
within the integration range). Note that the accuracy of 
the three-eighths rule is slightly inferior to that of Simp- 
son’s rule when the same tabular values are used; it is, 
however, useful to complete an integration over an odd 
number of sub-intervals, when the three-eighths rule can 
be used for either the first three or the last three intervals 
and Simpson’s rule for the remainder of the range. 

The above integration formulae use tabular values lying 
within and at the ends of the range of integration, and they 
are sometimes known as closed-type integration formulae. 
We can also develop open-type integration formulae, which 
involve only tabular values within the integration range. 
For example, if we use a parabola fitted to the tabular 
values fi, fe, and fg to derive an integration formula over 
the more extensive range (xo, xo +4h) we have* 


® Note that §=0 in the interpolation quadratic corresponds with the 
tabular point x, =%) +A. 


4.3 QUADRATURE FROM A TABLE OF VALUES 
zot4h 3 
| ᾿Ἰα)άχε [ fit 646: 40(0—1)4%f}d0 
ῳ Zo = 
= FM fifo 2) 


The correction term is + 14h5f((£)/45. Closed-type inte- 
gration formulae will obviously be more accurate than 
open-type formulae over the same range, as they make use 
of more information; hence closed-type formulae are usually 
used for quadrature, but open-type formulae are useful in 
solving differential equations. 


4.3. Quadrature from a table of values 

Formulae for quadrature from a table of finite differences 
can be obtained by integrating the interpolation formulae. 
Thus from Everett’s formula we have 


[Λα = 5p (fo-+fi)— 75 (80+ 4) + 


+ ὩΣ δίῃ + δι) ote saan δῦ if + δὺ fi} 4) O( h?), 


For integration over an extensive range we can add formu- 
lae of this type to span the range (using the results 
Sfo+ δῇ = df a of “} etc.), 


oe ‘Sleds sf tfitfot...+fnat afe & 


ee Paykel 
τη μδύη μδ[0) +55 (Hoya pd%fo) — 
<<.) Soe. OAT poor ΕΝ , gustan 
— spag0 δ᾽» — HB) + χε κερρί, δ» — δ} + OC), 
These formulae represent improved forms of the trapezoidal 
rule with corrections which involve the differences near 
ends of the range of integration; they may also be regarded 
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as giving the error of the trapezoidal rule in terms of finite 
differences of the integrand. 
Similarly, from integration of Stirling’s formula, 


a Re Pe Pe Scar’, 1 ἐν {19 
[2 J (x)dx==2h fit, +28f\ — 180. Ὁ + en" τοί 
= if +4fi+fe)— sph tae »τρδῦ τ Ο(ν3), 


giving the corrected form of at rule in finite dif- 
ferences. | 

The choice of truncation point in each of the above 
formulae must be made separately in the context of each 
new problem, but the mode of writing the formulae 
adopted here indicates the dependence on working interval 
of the remainder in case of truncation after sixth or seventh 
differences. In practice it may often be sufficient to work 
with only two or three terms. We often omit remainder 
terms when writing formulae; for in fact we seldom bother 
to evaluate them but rather infer their magnitudes from the 
size and sequential behaviour of the terms which we do 
take into account. Thus we judge directly from the table 
of differences where a formula ¢an be truncated without 
serious truncation error, and in so doing we select the 
degree of the approximating polynomial implicitly. _ 

The formulae of this section (used in truncated form 
with one or more difference terms) take into account some 
tabular values outside the range of integration as well as 
those within and so correspond with the use of higher 
degree matching polynomials; they are known as partial- 
range integration formulae, and for a specified working 
interval they will be considerably more accurate provided 
that f(x) is well-represented (and this can be judged from 
the size of higher differences). 


5 
Example 10. Evaluate | e*dx from the table 
0 
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oO Ss 1-0 1-5 2-0 2:5 
e* 1 16487 27183 44817 73801 12-1825 
x 80 3-5 40 4-5 5-0 


e~ 20°0855 331154 545982  90°0171 148-4132 

The integrand e* is a smooth function and should be easy to 
approximate except perhaps at the upper end of the range, 
where it is increasing more rapidly, The correct value of the 
integral to 4D is e&—1=147-4132. The trapezoidal rule with 
interval h= 0-5 gives 


5 
| e*dx=>'5(4 x 14+ 1°6487+ 2°7183+ ... +90°0171+ 
0 


+4x 148-4132) 
= 150-4716, with error 3-0584. Simpson’s rule with the same 


5 
working interval gives | e%dx<=(-5/3)(1+ 6°5948+ 5°4366+ | 
0 


+ 17-9268 + 14:7782 + 487300 + 40°1710 + 132-4616 
+ 109-1964 + 360-0684 + 148-4132) = 147-4628, with error 
‘0496. The corrected trapezoidal rule needs finite differences 
calculated near the ends of the range, given in the following 
partial table: 


—-5 -6065 «1549 0392 0113 
3935 +1003 0262 0043 

Ὁ 1-0000 «$52 0654 0156 
6487 1657 0418 ὉΙ2 

‘51-6487 -4209 1072 0283 

4:5 90-0171 22-9772 5-8643 1-4995 


ΠῚ 58.301] ΌῈ᾽Ο 49054 3-8059 9669 
5-0 148-4132 37-8826 6102 2 

2 6-2723 1-6022 
4-0686 


96-2787 5156 

5.5 244-6919 62-4582 15-9425 
5 

Hence | e%dx <= 150-4716—(-5/12)(77-3374— +5211) + 
0 


+ tS x 11/720)(19°7405 — -1330) — (δ x 191/60480) x 
x (5°0391—-0340) + (-5 x 2497/3628800)(1-2846 — -0085) 
= 147-4132. 
Note: (i) that the simple trapezoidal rule is most seriously 
in error near the upper end of the range where εὖ varies rapidly 
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. 5 
with x and there is a contribution 1-9464 from | e*dx to the 
4 


whole error 3-0584; (ii) that even with so large a sub-interval 
as h=(-5 the parabola fitting of Simpson’s rule gives very rea- 
sonable accuracy when the function is smooth in this sub- 
interval; (iii) that a standard method of checking an integral 
evaluation is to repeat the calculation using a smaller interval 
(say half of that used initially) and when the two results agree 
to the accuracy required we can be confident both that the 
original interval was suitable and that the calculation is free 
from blunders; (iv) that when greater accuracy is needed the 
finite difference formula should be used; (v) that when the 
behaviour of the integrand varies appreciably from part to part 
of the integration range the most effective procedure may be to 
use working intervals of different length in different parts of the 
range so that the local accuracy of working is roughly uniform; 
(vi) that blunders are specially likely near points at which the 
working interval is changed, because of the interruption of a 
smooth routine, and as a safeguard the two calculations should 
be overlapped to give a check. 


1 
Problem 4. Evaluate | tan νὰ dx to four decimals. 
0 


Our first step might be to tabulate the integrand and its 
differences at interval h= 0-1 (see next page). If we evaluate the 
integral by a normal application of Simpson’s rule with work- 
ing interval k= 01 we have from the table, 


| jan Vide (0+ 4X 3272+2x -4796+ ... +4 139454 


+ 1:5574)=0-8536. 
This seems straightforward enough, and as a check on our 
working we repeat the calculation with h=0-05 and we find 
the value 08553. This poor agreement is not due to a blunder. 
It can be seen that differences near the top of the table are large 
even at high order, and this is related to the fact that the deriva- 
tives of tan νι are all infinite at x=(. In fact all integrands 
which behave like x’* or x-“ near a terminal point x=0 of an 
integration range (or like |x—a|+t near x=@) are troublesome 
7 
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x ὟΣ] tan ν' ἃ ὃ δὲ §3 δὲ 
0 0 0 | 
ahs 1748 
. *3162 *B272 —174 
: : 1524 1528 voi 
"4472 -4796 — 220 —13 
: 1304 - 145 δι 
.8477 ‘6100 —75 wn 
: 1229 7 62 Ἢ 
. "6325 ‘7329 --Ἰ -- 
: 1216 - 88 P 
“ἢ 07] “8545 + | -- 
. 1241 30 
6 “7746 “9786 55 —5 
1296 ms 25 Σὰ 
. ‘8367 110 82 { 
1376 31 
"8 ‘8944 1:3458 111 0 
1487 31 
“9 ‘9487 1°3945 142 
1629 


1:0 10000 1°5574 

even though their integrals converge. This kind of difficulty can 
be sidest d in various ways. 

εἰς We cai use an open-range formula for one set of intervals 
adjacent to x=( (allowing us to ‘stand back’ a little from the 
trouble level x=() and Simpson’s rule in the remainder of the 
range. Thus with A=0-1, 


1 4 1 
[tan vede= | tan ya de+ [tan vids 
0 0 . 


~(4x +1 x $)(2 x *3272—+4796+ 2x *6100)+ 
+(-1x 4)(-7329+ 4 x 8545+ 2x *9786+ ...+ 
+ 1:5574) 

= 08582. 

5 correction terms — f(#)( £) x 10-5/90 for each application of 
Εἰκεδηρωςς rule can be Hic negligibly small, but there is 
doubt about the term + 14f(( £) x 10-*/45 for the initial open- 
type formula. Probably the — way of resolving this doubt, 


———VOOoe a 
ea = 
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and at the same time checking the calculation, is again to repeat 
the process with working interval 0-(05, and the ‘improved’ result 
is 0-8569. This decrease in working interval has clearly de- 
creased the error from the open-range part of the integration 
but we cannot be sure that it has eliminated it. All that we can 
reasonably assert at this stage is that the integral lies between 
0°8553 and 0-8569, as the correction terms are of opposite sign. 

(b) In the neighbourhood of x=0 the integrand behaves as 
/x, and hence if we compute the integral of (tan /x—4/%) 


and subsequently use the result | νὰ dx= 2/3, the numerical 
0 


work should be a good deal easier. Using values from the table 
above, we have from Simpson’s rule with :=0-1, 


1 1 
(tan /%— ν᾿ ἀ)άχεεξ" 1896, so that | tan ναὶ dx=>-8562. 
0 0 


When the integral is re-evaluated using working interval 
h=0-05 the result is again 0-8562, and this may be taken as 
correct to 4D. 

(c) An important technique for dealing with awkward 
integrals is to convert them into differential equations in which 
the dependent variable is free from the singularity that has been 


= 
causing the trouble. If we take | tan /u du=g(x)y(x), where 
0 


(x) is a function to be determined and g(x) is selected as having 
a singularity related to that of the integrand, we obtain a dif- 
ferential equation on differentiating. With g(x)=/x (same 
singularity as tan ναὶ at x=()) we find Qxy’+ y=2,/% tan ν 2; 
with g(x)=° (when dg/dx has same singularity as integrand) 
we have 2xy'+3y=2tan νι να. For the solution of these 
differential equations see Numerical Solution of Equations. 

(d) A change of variable often serves to reduce a difficult 
integral to more tractable form, In tie case, if we take /%¥=u, 


1 
tan νὰ dx=2) utanu du, 
0 } 0 

and numerical evaluation is considerably easier. 


Note: (i) that for integration it is normally sufficient to 
work with four decimals to produce a result correct to 4D; 
(ii) that it is worth tabulating the differences of tan +/% as a 
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check on the calculated tabular values; (iii) that care is needed 
in working with differences where there is a sign change that 
might be overlooked ; (iv) that the variations in fourth differences 
are not so large as to indicate an error in calculation; (v) that 
a difficulty involved in the pursuit of method (a) to smaller 
working interval for the open-range section is that f@(£) will 
at the same time increase rapidly; (vi) the methods of (b) and 
(c) are sometimes known as ‘subtracting out the singularity’ and 
dividing out the singularity’, respectively. 


oa 
Problem 5. Evaluate z2to 6D from the series 72/6= > (lm). 
m=1 


The process of summing terms of a series is related to that of 
integration. In general we need distinct formulae for the two 
types of process, but when the corrected form of the trapezoidal 
rule is expressed in terms of derivatives of the integrand (in 
place of finite differences) we obtain a formula serving both 
processes, and this is the Euler-Maclaurin formula*, 


Zotnh 1 a adh 
κ᾿ fls)ds=h{ Sor fit fot oe t faut stn) — RU, + 
+ gle —I8)— aa —19)+ 
PHM ΠΡ). vee 


For the summation of terms of a series the Euler-Maclaurin 
formula can be rewritten: 


fotfit ... ἘΔ 


A . 
mje) dt ot + GML) Us fe) 


ΝΡ 1 
Ἐιβσθ- BS) πο 08"- FP 0 + 


® For a derivation of the Euler-Maclaurin formula see for example, E. 
T. Whittaker and G. N. Watson, A course of modern analysis (Cambridge 
University Press), chapter 7, 83; it can also be derived from the trapezoidal 
pr schun for an extended range by substituting for finite differences in terms 
of derivatives. 
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Hence, when the general term of a series can be written f(m), 
where f(x) is a function possessing the necessary derivatives 
and which can be integrated, we can form sums of terms of the 
series. In the present case f(m)=m-®; thus f(x)=a-* and 
f™(x)= ss 1)*(k+1)!/x*+2, To sum the series take h=1 and 
x=1, when 


Th 
ee ee οὐ 
aa Ι ' τὰν - it ἢ te ety) 


tif. 98 4 


and if we take the limit of the sum as πω, 


co 
SL a SS A Sl 
Dam it ant Est DE "8 


These terms decrease very slowly, and we should have to 
accumulate far too many for the required accuracy. This is 
because f(*)(1)=(—1)*(k+1)!, and it is apparent that a much 
more effective process can be obtained if we sum the first few 
terms directly and apply the Euler-Maclaurin formula to the 
remaining terms. ‘Thus 


ps a] 
20h) τ SST ee ee 
Pp τοῦ 5105 6.108 80.105 1 45.107 ᾿" 
= +40000000-+ -00500000-+ 00016667—-00000033-+ “00000000 


= 10516634. 
By direct addition 


>. 4 
>. 153976773 
1 


and hence 72=6x 1°64493407=9°869604. The obvious check 
is by repeating the process, taking a different number of initial 
terms for direct addition, and this must be repeated until two 
successive answers agree. 

Note: () that this form of the Euler-Maclaurin formula does 
not imply that it is an infinite series, and for many integrands 
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it will prove to be an asymptotic approximation in which the 
accuracy improves with each successive term only up to a cer- 
tain stage and thereafter decreases progressively; (ii) that the 
Euler-Maclaurin formula can be used in its first form for inte- 
esa of a specified function provided that derivatives can be 
ormed and evaluated with reasonable ease; (iti) that we must 
retain eight decimal digits if after adding 9+ 4 terms (each with 
possible round-off of “δ in the last digit) and multiplying by 
6 we are to achieve a result correct to 6D; (iv) that to obtain 
6D accuracy in πῇ by direct addition of terms we should have 
to add some twelve million terms of the series! 


4.4. Double and repeated integration 
Double (and other multiple) integrals can be evaluated 


by use of the foregoing general techniques applied to a 
repeated integral. Thus 


[i [ope y)dxdy= |" {{{ Ste. y)ds\dy, 


h 
and if we take F(y) =| fm y)dx, we have from Simpson’s 
rule τ᾿ 


[᾿,{1|ργακάν 
<7! F(—1) + 4F(0) + F(h) 
-" ΓΌΩΝ | * fs Ode | fl iyas| 
yh sal a+4fo-1+fi-] + FhLf-10+ 4foo + fio] + 
+3h [f+ 4for +ful} 
= iiffint faut faatfiat sfiot Mort 4f-s0+ Moat 


16foc}, 
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where fy= (th, jh). In principle this is a perfectly satis- 
orm of Simpson’s rule for a double integral, in 
which the value of the integral is given in terms of tabular 


factory 


values at the junctions of a square mesh of side-length 2h, 
as indicated in figure 4. However, it will obviously be 
laborious to apply, and we should prefer a formula with 
fewer terms on the right-hand side if this could be obtained 
without undue loss of accuracy. 


ri 


If 1-1 . ΠΥ 


ἘΞΕ τοὸοῦὸῊ eS ὑπ ὑπ πῖστα — SS ὅσσ ὅσ, ὑἱἱἱ an anal 


fir 


Fic. 4. 

One simpler formula, in which the x, symmetry i 
preserved but from which the relatively distant πολιῶν 
ae Su, f-u, f--1, fi-1 are omitted, is of the general 
orm 


h fh 
| [i fe» noddy=afto+ i fio+ far +f-s0+fo-s, 


where a and 6 are constants to be determined. These con- 
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stants can be found easily if we make use of the idea of 
polynomial approximation in a slightly different way, and 
in particular of quadratic approximation which we have 
already found to be effective in integration. If the formula 
that we are constructing is to correspond with quadratic 
matching of the integrand f(x, y) then it will integrate 
correctly (i.e. without any error) all quadratic functions 
of x and y, including of course linear functions and con- 
stants. Moreover, since integration is a linear process, we 
can be sure that all quadratic functions will be integrated 
correctly if this is true independently for «?,x (and 
hence by symmetry for y?, y), xy and 1. 


When f(x, y)=1 we have from the formula 
hh 
| [ 1ἀχάγτεκα- 4; 
J-al - 


this double integral can be evaluated directly and has the 
value 4h?, and hence the formula integrates constants 
correctly if a+4b=4h*. For f(x, y)=x, 


(rh Ah 
o=| | xdudyax0+bx0, 


so that the formula integrates x, and indeed xy and all odd 
powers of x or y (and products odd in x or y) correctly. 
For f(x, y)=2*, 


h th 
(2h8/3)2h= | J rdxdya x 0+5 x 2h, 
and the formula will integrate x* correctly (and also y?, by 
symmetry) if b=2h?/3. Hence the formula will integrate all 
quadratic functions correctly if a= 4h? -—4b=4h?/3, and we 
have the double integration formula based on quadratic 
matching, 
h ph vee 
PP feo ardedy Fifa fio + far+f-10+ fos) 


= 4(area of element)(2fo0 +f10 +fo1 +f-10 +fo-1). 
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Other formulae of this general type can be obtained in a 
similar way, and these can be tailored to fit any special 
requirements of the problem in hand; indeed, the deriva- 
tion above has been given in detail to demonstrate the con- 
struction of a formula as well as for its own importance. In 
each case we can calculate a correction term provided that 
a Taylor series (in x and y) is available, and the correction 
terms for the formulae of this section are 


ie —(h8/45){ fexral ξ, ἢ) t+fuvyv( é, »)} 
— (δ 45){{π::«( ξ, ἡ) -- 5ϑσεννί ξ, 1) ἘΪννννί ξ, a}, 


respectively, where (€, 7) is a point of the square domain 
of integration and fzrzzz=C4f/ox', εἴς, 


4.5. Gaussian quadrature 

Our previous integration formulae have all been based 
on the use of equally spaced tabular values of the integrand. 
Such formulae are convenient for integration from con- 
stant interval tables, and also for integrands of given func- 
tional form (as the tabular points can be chosen for easy 
calculation), but they are not the most accurate type of 
formula obtainable. If instead we specify in advance the 
number of tabular values that we wish to incorporate in an 
integration formula but not their position, we can construct 
a formula by finding the tabular positions within the range 
of integration for the most accurate integration formula 
obtainable with the specified number of values of the inte- 
grand. ‘This approach is due to Gauss. 

Thus, the most general form of two-point formula of 
this type is 


[7 florax=nfafiah) + 54080), 


where the weighting constants a and b, and the position 
factors « and f are constants to be determined. The approxi- 
mation symbol = has been used in anticipation of the fact 


4.5 GAUSSIAN QUADRATURE 


that we shall now obtain formulae based on polynomial 
approximation (using the device of the previous section). 
The simplest formula of this type* is symmetrical about 
x=0 and uses values of the integrand calculated at points 
--ah; for this case B=—«. We now obtain the values of 
a, b, and « by stipulating that the formula will integrate 
1, x, and x? (and hence all quadratic functions) without 
error. 

The constant 1 will be integrated correctly if 2h=h(a+ δ); 
and x correctly if 0=h(axh—bah); and x* correctly if 
2h3/3 = h(ao®h? + bo?h?). Hence a=b=1, and o2=1/3; and 
we have the Gauss two-point integration formula 


[7 flrax=n(f0h/ v3) ἘΚ —H1V3)} 


this formula (like Simpson’s rule) integrates powers up to 
x3 (and all odd powers of x) correctly. 

The Gauss two-point formula integrates quadratics 
correctly and corresponds to quadratic representation of 
the integrand f(x). The error may be calculated using 
Taylor expansions, and the leading correction term is found 
to be +(h5/135)f( €), which is very similar in magnitude 
to that for Simpson’s rule but is of opposite sign. It should 
be noted that the Gauss two-point formula achieves as 
good accuracy as Simpson’s rule in spite of the fact that it 
employs only two instead of three tabular values; this is 
the advantage of Gauss type formulae, but it must be com- 
pared with the disadvantage that these values must be 
determined at awkward points. Thus when we are working 
by hand the inconvenience of evaluating the integrand at 
the Gauss points will more than offset the improved 
accuracy in general and Simpson’s rule (say) will be more 
effective; but when we are working with a high-speed 
computing machine which can handle any tabular points 


® For special applications we could easily develop non-symmetrical 
formulae of this type. 79 
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with equal ease the Gauss formulae will have real ad- 
vantages. 


The Gauss three-point integration formula 


[i fone 
= 5h(5f (-1,/ 3 + 8.00) Ἐ 57 (h | 3)} + (h7/15750)f £) 


can be obtained in a similar way, as the formula of type 
h{af( — ah) + af(ah) + bf(0)} corresponding to cubic approxi- 
mation to f(x). 


EXERCISES ON CHAPTER FOUR 


1. Derive the half-way differentiation formula 
Lf -Δ. _ 1 99 3 ΜΒ os 
= 8, χη, 540° a 716 g fy, HOM). 


Find f'(:55) from the following table for f(x)=e~*, and check the result 
by interpolation. 


x 3 “4 5 6 7 ‘8 
Six) 74082 67032 -60653  -54881 -49659 -44933 
2. Obtain the differentiation formulae, based on tabular values, 
2hfo' =fi—S-1— bP" f(), ' 
12hfy' = Sat Bhi Bf that aH FO). 
3. (i) Derive the three-eighths rule for integration 
ah 
| Sx)dx = ἐπ 0 30. +3fa+Ss) 
J 0 
using the matching polynomial P,(x), and show that the correction term 
can be written as — sng (£). 
(ii) Derive Weddle’s rule 
ἢ 
[eres a aahlfo-t Sh that Ofathetist+he) 
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EXERCISES 
from the polynomial P,(x). (To obtain this form a small term must be 
neglected.) Find the correction term. 
(iii) Compare the accuracy of Dufton’s rule 
104 
5 
| f(x)dx τος Sif +Sa t+ Sethe) 
o 
with that of the trapezoidal rule and Simpson’s rule. 
4. Derive the integration formula for use with divided differences 


[ω] 
| aia = (x, — Xp) f(%o) +44 —%0)* Loe) -οξίαι —%)*Lo% 2], and 


show that the correction term is [Ὁ --αρ) αν —x,)(% πολ  ΓΧΧ λα ]ἶχ, 
To 


1-5 
J(x)dx from the following table for the error function 
0 


τιν πεῖ 2, comparing the results obtained using the trapezoidal 
rule, Simpson's rule, the three-eighths rule and Weddle’s rule. 


5. Evaluate 


x Fix x F(x) x 2) 
000 Aste “625 65632 1-250 JB. 
125 -79168 ‘750 60227 375 31002 
‘250 -77334 “875 54411 1-500 25903 
375 74371 1-000 -48394 
500 70.413 1-125 42375 


6. Integrate 1/(10-+x*-++,y*) correct to 5D over the square region with 
corners (+1, +1). [Carry out the integration with the region divided into 
one, four, and then nine elementary squares, and compare the values ob- 
tained from these three different working subdivisions of the region.] 


7. Derive the repeated integration formula 
h ἃ; 

1 1 i 7 1. | 
[ as S(Odt = 5h (1-+3p +75°° fo— {apt fo — gag? fot .) 
0 0 
Hence deduce the repeated integration formula 


an τ Te 31 | 
| Wicca (Caprice 540° ον Baap fot -) . 
- Γ 
8. Construct a double integration formula of type 
h h 
[ dx| f(x,y) dy ea foot δ ti Hoa th). 
-ὰ 


J-h 


What kind of polynomial matching does it correspond to? Give a Taylor 
series estimate of the correction term. 
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CHAPTER FIVE 
Method of Least Squares 


Finite difference processes may be used effectively with 
tables of well-represented functions because the corre- 
sponding differences decrease rapidly with increasing 
order. Thus higher differences of mathematical tables 
(with reasonably short interval) tend to be small and regu- 
lar in behaviour; and indeed we detect occasional blunders 
by the relatively large disturbances they cause in the other- 
wise systematic behaviour of the difference table. The same 
processes will generally be ill-suited for use with tables of 
experimental results, as these will include random fluctua- 
tions that often dominate the higher differences. Most 
experimental values will be in error (by an amount which 
exceeds normal round-off), and the effects of these errors 
in the higher differences will be so intertwined that there 
is no hope of detecting (let alone of correcting) them by the 
method of §2.5. We must develop some other way of 
᾿ς handling experimental results. 

So far we have selected our approximating polynomials 
to agree exactly with the specified tabular values over a 
chosen range of tabular points; thus we tolerate some error 
of approximation between tabular points, though none at 
the points themselves, But when the tabular values include 
random errors there is no advantage to be gained from fit- 
ting them precisely, and we would prefer approximating 
polynomials that follow the general trend of a table without 
reproducing all the local fluctuations. This relaxation of the 
accuracy with which we fit individual tabular values in- 
volves a reduction in constraint on the approximating poly- 
nomial and a consequent increase of range over which a 
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5-1 FITTING STRAIGHT LINES TO DATA 


polynomial of given degree may be used in approximation. 
Matching processes which involve some averaging of small- 
scale local fluctuations are known as curve-fitting. 


5.1. Fitting straight lines to data 


It is often obvious from graphical inspection that the 
measured values, say γεξξ γίχι), obtained from an experi- 
ment lie roughly along a line. These results apparently 
define a linear relationship of the form 


y = a t+ax=P,(x), 
where the linear representation P;(x) is not a Newton poly- 
nomial (which would have zero error of approximation at 
tabular points) but is some sort of averaged representation; 
we want a procedure for selecting such a ‘line of best fit’ 
or ‘average line’ through the points. It should be noted that 
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there is no ultimate best line, unless the points are collinear, 
and we shall obtain different ‘best’ lines according to the 
criteria of selection we adopt. A good line will pass close to 
most points, and we can select such a line if we adjust ao 
and a; until the sum of the squares of the residual distances* 
¥i— (a+ a1%;) is a minimum. The use of (distance)? rather 

simple distance provides a weighting towards distant 
points so that points already near the line exert relatively 
less influence and contributions from points on opposite 
sides of the line reinforce. This criterion of least squares will 
ensure that the distribution of points as a whole lies closely 
about the best line; it is not the only criterion that could 
be used, but it is convenient in many applications. The 
method of least squares can best be illustrated by an example 
of linear approximation. 


Example 11. Use the method of least squares to fit a straight 
line to the following data. 
x 0 1 2 3 4 5 
y O82 346 5692 857 11:07 13-64 
The sum of squares S of the residual distances relative to the 
line y=ao+ a,x is 
S= (0-82 —a)?+ (3-46—ao— ay)? +-(5-92—ap—2a)2+ ... | 
+(13-64—a9— 5a). 
In order that S may be a minimum for variations of a and a; 
(1.6. for variations of the line) it is necessary that 
os 


o=— 


a 
= — 2{(0-82—a0)+(3:46—ap—ay)+ ... +(13°64—a—5ai)}, 
0-5 
= — 2{(3°46—a9—a,)+ 2(5-92—ap—2a))+ ... 


+ 5(13°64—ao— 5ay)}; 
and hence that af } 


* Note that these residual distances are measured parallel to the y-axis 
or ease, but are proportional to distances from a line 
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θαυ" 15a; = 43°48 
1549+ 55a,= 15349. 
Thus the line of best fit in the least squares sense is 
y == 0°85+ 2-56-. 
te: (i) that we have not proved that this line corresponds 

to pects of S, though this is not difficult to do; (ii) that the 
approximate values at tabular points are 0-85, 3-41, 5°97, 8°53, 
11:09, 13°65; (iii) that a measure of the effectiveness of the 
approximation is provided by the root mean square of the 
remaining residuals | 

{ξ [(-03)2+ (-05)?+ (-05)?+ (-04)?-+ (02)? + (-01)?}}#= 0-037; 
(iv) that we cannot strictly determine the accuracy of 
approximation as we do not know the correct tabular values; 
(v) that we can fit a Newton polynomial Pj(x) to just two tabular 
points, but that a least squares polynomial P(x) can be fitted to 
any number of points. 
5,2. Polynomial approximation by least squares 

In very much the same way we can find er 
approximations P,(x) of any degree for given data; t 1e 
hain of degree is in no way regulated by the number of 
data points, though the accuracy of approximation will 
depend on our choice. The procedure can conveniently be 
described in terms of the algebraically simple case of fitting 
a parabola P2(x) = ao + a,x + aex® to the n+1 tabular values 
f(x) (n>2, and ¢=0, 1, ... 2). | 

We must choose the coefficient parameters ao, ai, and 
dg so as to obtain a minimum value for the sum of residualr 


s-> £ flees) — ao — ares — ar9}?, 


A set of necessary conditions for S to be a minimum is 
that 0S/éa;=0 for k=0, 1, 2; and hence 


=2> {f(%)-ao—aix—anxaxt=0 k= 0, 1, 2, 
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nn n 
If we take 5.5 > xf and v,= > xf, these equations 
i=0 i~0 


can be written 
Sod + 5141 + S2a2 = U9 
δι + 5561 + 5342 = Ὁ] 
52d + 5301 + S42 = Va, 


and are known as the normal equations of the system. It can 
be shown* that the normal equations always have a unique 
solution for ao, a1, @2 and that this solution does correspond 
to a minimum of S. 

The derivation of the normal equations should be laid 
out systematically, as in the following table: 


1 Ὁ 

1 Δ < oct 1 fi 1 δὴ ἢ wif 
Lo * κα of fe κε καἶ 
χᾷ 


Ὁ 81 52 88 54 vo YU ὍΣ 
The values x; and f; are given, and may be entered forth- 
with; their powers and products are then calculated and 
entered, and the columns are added to obtain s and v sums, 
Example 12. Fit a least squares polynomial to the following 
experimental data (which may be in error to the extent of about 
ὍΣ per value): 
” 0 1 2 3 4 5 6 
y ““Ὃ 37 ‘71, 98 1288 41658 1-69 
A substantial saving in labour can be obtained from the 
change of origin ¥=*— 3: in terms of % the tabular points are 


® See Milne, Numerical seas gece 1949). 


5-2 APPROXIMATION BY LEAST SQUARES 
symmetrically arranged about #=0, and if the tabular interval 
is uniform the odd sums are zero, i.e. 1=0, 3= 0, we 

We obtain the normal equations for a /inear approximation as 
follows: 


1 Ν: ee: ἂν 
1 —3 9 -- Ὁ] "03 
1 —2 4 ‘37 - 74 
1 --1 1 71 - 11 
1 0 0 98 0 
1 1 1 1:28 1:28 
1 2 4 1°52 3°04 
1 3 9 1 69 507 
Υ 0 828 6°54 7°97 
Hence Ταρ- θαι] =6°54 
OQag+ 28a,;= 797; 
and y = 0°934+ 0°285* . 


A measure of the effectiveness of this approximation is provided 
by the r.m.s. (root mean square) of the remaining residuals ; the 
remaining rms residual is 0-06, and as this is substantially 
larger than the anticipated experimental error it appears that 
this linear approximation is inadequate. _ er 

The normal equations for the quadratic approximation fol- 
low as: 


a2 | μὰ μὰ μεὶ μαὶ μὰ μεὶ μεὶ 


28 0 196 654 797 2467 


7ao+ θαι + 28a2= 6°54 
Qao+ 28a; + 0ag= 7:97 
2849+ 0a) + 196a2= 24°67. 
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Hence y = 1005+ 0-285 ¥— 0-01 23 
= —0-012+ 0°393 x—0-01g x*. 
The rms residual error is now 0-01, which is well within the 
range of possible experimental error and must be regarded as 
satisfactory. 


Note; (i) that we now rely on our judgment in choosing the 
degree of approximating Sclencuntsl that is appropriate, and 
there is nothing to prevent us from finding a ‘best straight line’ 
through points obviously ill-suited for linear approximation; 
(ii) that the root mean square residual is a valuable way of 
assessing an approximation by least squares; (iii) that there is 
an effective reduction by one in the number of normal equations 
to be solved when the tabular points are arranged symmetric- 
ally; (iv) that the third decimal digits are needed as subscripts 
(but cannot genuinely be written as ‘correct digits’) in order to 
preserve accuracy. 


Experimental results are often plotted on log-log or 
log-plain graph paper to bring out relationships of power 
law and exponential type. The same procedure can be 
carried out much more accurately using least-square linear 
approximation. In the case y=ax", we have log y=log a 
+n log x, and if we put X=log x, Y=log y, and A=log a, 


we have, 
Y=A+nX; 


A and n can be determined now by a least-squares process 
applied to the values logy; at points log χε. Similarly, if 
y=ae*, wetake Y=log y and A=log ato give Y= A + dbx, 
and a least-squares process can be applied to the values 
log γι at points χε. 

It should be emphasized that the least-squares approach 
is bound to yield an approximation,* though we shall get 
a useless approximation if we select the degree of the ap- 
proximating polynomial foolishly. It follows that we must 
always test the significance of an approximation by a final 
calculation of the rms residual, or in some other way. 


® Compare the comments on — theories of 83,7 (iii). 
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In some applications the tabular values may be regarded 
as of unequal importance or weight; in such cases a weight 
Junction w(x;) can be associated with each tabular point, and 
the least-squares process will then be carried out on the 
weighted sum 


S= > e( 4) {f( 1) — Pn(xi)}?. 


Problem 6. Find a least squares quadratic approximation for 
the continuous function e? over the range (— 1, 1). 

We might select an arbitrary set of tabular points spanning 
the range (—1, 1) and apply the methods of this section to the 
corresponding set of values for e%. But it will be obviously 
superior to use the known values of e* right across the range by 
using integration in place of the summation process necessary 
for sets of discrete values. Thus, in place of the sum of squares 
of the residuals we take the integral 


1 
S= | or (49+ ayx+ agx*)}2dx, 
which is a function of ap, a; and ag. A necessary set of conditions 


for a minimum of the function S(ao, a, ag) is that dS/dap=0 
for k= (0, 1, 2: hence 


ib. hill f {e%—(ag+ ayx+ apx")}dx 
Cao 3 1 ι 
0= τ - {e%—(ao+ a1x+ agx*)}x dx 
as gas κοὐ, 22 
= -9 | x, (@o+ ayx-+ apx*)}x? dx, 


After integration, we obtain the set of normal equations 

3ag+ az = 3:53 

a;= 1-10 

5a9+ 3a2= 6°59, 

and the quadratic approximation 
y = 100+ 110+ 0-53x%. 
Note: (i) that the numerical value for e has been substituted 
immediately as we are concerned with numerical solutions; 
»" | 
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(ii) that the error does not exceed about -03 units except at the 
very ends of the range, and that the rms residual error is -03; 
(iii) that the approximation differs from the first three terms of 
the Taylor expansion near x=(, but is more accurate over the 
specified range as a whole; (iv) that this process of matching a 
continuously defined function is known as continuous approxt- 
mation by least squares. 


5,3, Data smoothing 


Raw experimental results provide our only source of 
contact with the physical world; we should on no account 
‘adjust’ these results capriciously, either to improve an 
‘apparent pattern of dependence’ (which may not be a 
genuine one, in any case) or to provide a better level of 
agreement with a favoured theory. Experiments which 
yield consistently poor results may need refining, or they 
may be based on an inadequate conception of the pheno- 
mena involved; in either case we can improve the apparent 
quality of the results by smoothing processes, but we will not 
improve their significance (though we may appear to do so!). 
Thus data smoothing is not an alternative to good experi- 
mental technique. However, it must be appreciated that 
even good experiments yield tables of results with enough 
experimental scatter to produce sizable spurious fluctu- 
ations in the higher differences. Some initial smoothing 
may be necessary before such tables are used in numerical 
calculations, especially for processes like differentiation 
which are very sensitive to local error. 

In smoothing a table of values we do not attempt to find 
a single approximating polynomial for the whole table; 
instead we use repeated local matching with low-degree 
polynomials, and each polynomial is used only to adjust a 
single central value (except at the ends of the range). Thus 
we assume that the table is Jocally suitable for approxima- 
tion by polynomials of low degree, apart from the experi- 
mental scatter, but we make no assumption about large- 
scale approximation of the νὰ 
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5.3 DATA SMOOTHING 

A wide range of formulae is available for smoothing data 
at uniform interval; we shall derive a quadratic five-point 
formula. We shall suppose that a run of five values is suf- 
ficient to correct for random errors in the measured value 
of the central member, and we shall take the tabular 
interval as the unit of length and the central point as 
origin. The normal equations for fitting a parabola 
ao + oP aox® to the values f-», f-1, fo, fi, fe by least-squares 
are then 


2 
5ao+ 10a2=> fi 
2 
Ξ 
10a, => xfs 
-2 


2 
10a9 + 34a2= ΣΝ. 
—2 


For use in the interior of the table we want only a corrected 
value for the central point, that is the value ao; and we 
have on solution 


35a0= 35fo —3(f-2—4fa+t Sfo—4fitfr) 


= 35fo— 35%. 


Hence the corrected value* at x=0 is yo = fo—38%f0/35. 
For use near the head or foot of a table 


1 1 
y-2=f-2- apo fos 42 == fo— ao so 


2 2 
yafatzdYo, γι =fit x25'fo- 


In selecting the degree ἡ for the smoothing process we 
should take small so that individual calculations are easy, 
and yet we must so judge our choice that differences of the 

® For many purposes it will be sufficient to use the approximate form 
Vo ee Soe  ἢ.12.Ψ 
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true function of order greater than m are small. In fact the 
choice may often be made from the fact that the first π 
differences are reasonably regular, but the (π 1.18 dif- 
ferences are irregular and have mean value near zero. The 
amount of smoothing increases with the number of values 
fitted, and decreases with the degree of the approximating 
polynomial. Obviously there is great scope for judgment on 
the part of the computer as to what is random and what 
genuine; sometimes the application of successive smooth- 
ing processes may be justified, but repeated smoothing will 
ultimately reduce any set of data to a straight line. 


EXERCISES ON CHAPTER FIVE 


Ped ny the method of least squares to find the best straight line fitting 
t ta 
so 3 2 3 4 5 6 7 891 
gy 15 14 13 13 12 10 6 8 5 4 2; 
2. Find the best curve of type vax" fitting the data 
« 013 033 108 3:65 4-98 
y 3 + 6 9 10, 


Note on Operational Methods 


Most readers will have some familiarity with operational 
methods through their use in solving linear differential equa- 
tions with constant coefficients. It is ible in a purely formal 
way to dissociate the or of differentiation (or integration) 
from any particular function f(x) and to establish a formal 
algebra for the operator D=d/dx, which is similar to (though not 
identical with) the algebra of real numbers. 

In precisely the same way we can regard the operations of 
differencing as distinct from their application to a particular 
function, in which case 4, 7 and ὃ may be regarded as differenc- 
ing [pace and μι as an averaging operator. A formal algebra 
can be established in terms of these and other operators, and 
this operational algebra can be used in deriving a very wide 
range of formulae for numerical processes. Operational methods 
are used in many books because they allow of such easy develop- 
ment of Sammie lac, but they have not been used here as they 
tend to leave the significance of the various steps obscure to the 
inexperienced reader. We have advanced more slowly, in the 
hope that the reader will carry a clear understanding of the 
reasons behind our development of processes and of our atti- 
tude towards numerical work; he can then take up the opera- 
tional approach with help from a more advanced text, 
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Answers and Notes on Exercises 


Chapter One 


1. (i) 3, 6; about 1 in 2x 10°, (ii) 30, 3; about 0-05 per cent. (iii) 0, 5; 
about 1 in 60,000. (iv) 4, 3; about 0-4 per cent, and note that the fi 
zero is needed to define this accuracy. 


2. 1:38 and 1-37: note that the rounding error can (just) exceed half a 
unit in the last digit when a number is rounded digit by digit. 


(i) 0-5; actually 0-47+7 x-005, and hence the maximum error in the 
second decimal from round-off is +34, so that this digit is uncertain 
(the probable error is considerably smaller, but exceeds -005). 
(ii) 83-5; maximum error is about (9-37+8-91) x -005=0-09, and the 
second decimal in 83-49 is unreliable and the error could affect the 
first decimal. (iii) 25-6; maximum rounding error is 0-048 after 
division and the second and third decimals are unreliable (clearly the 
answer cannot be as accurate as the denominator 0-279). (iv) 17-760; 
the decimal accuracy of these numbers varies and as the middle num- 
ber is known only within -0005 there is no justification for retaining 
more than three decimals, maximum error is about -0006. Note that 
maximum errors can be found using the binomial expansion, as 


(atey(@18)=$( 185) (118) "= {te(,+3) + 


4, (i) 0693144 with maximum rounding error of 1 in last digit, (ii) 
0693143 with maximum rounding error of 4 in last digit; this illus- 
trates the value of carrying an extra working digit for some calculations. 


5. (i) 1+”, for the range —1<x<-1; (ii) i+x+4x, for the range 
3 <x 3: (iii) 1+a+4x2+424, for the range --ϑ χα «5. 


3 


᾿ 


Chapter Two 
2. Compare with formula for differentiation of products, etc. 
3. Use induction. 


4. Third differences are constant: hence data fitted by a cubic which is 
obtained from an interpolation formula as 75 +521x/6—25x*/2 + 2x°/3. 


5. f(1-5)=1-142892. Transposition is a common blunder; the adjusted 
fourth differences are about as regular as can be expected, allowing 
for round-off. 
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Chapter Three 


1. 


3. 


4, 


Use Newton-Gregory interpolation with forward differences near 

head of table and backward differences near foot. Differencing of final 

ay will uncover any errors, so that there is no need to give solution 
ere. 


. (i) Suitable for linear interpolation as 6°f < 1 over whole table. (ii) Not 


suitable, as δὲ rises from 15 to 29. Subtabulation to half-interval will 
reduce second differences by factor of about 1/27; they are still too 
large to be neglected. Subtabulation to one-third interval will reduce 
second differences by factor of about 1/32, when they may be neglected 
over whole range (as being <4 units in last decimal). 


f(4-2)=0-238095. Throwback of fourth on to second differences is 
not good enough because of the large interval and number of decimal 
digits; use Everett to sixth differences. 


§¢f is about 10, and makes a contribution to Bessel’s formula of about 
one unit in the last decimal over much of the §-range! Values for the 
tangents using Bessel’s formula are 2-36158, 2-40827, 2-46888, and 
using mean differences are 2-36181, 2°40831, 2°46866, respectively, 
and thus errors of more than twenty units in the last decimal are possible 
near 0’ and 60’! In fact the values using mean differences for 67° 3’ 
and 67° 57’ are correct to 3D only!!! and that for 67° 27’ to 4D! It 
is hard to see why the author retains mean differences for 67° but 
rejects them for 68°; and the reader may find it instructive to investi- 
gate other short sets of tables and to compare them with Chambers's 
Shorter Six-figure Mathematical Tables. 


. First reaction is to use the Newton-Gregory forward difference inter- 


polation formula as we are working near the head of a table, but we 
find that the terms decrease too slowly for effective calculation. 
(N.B. if we take many terms this will involve many tabular values, 
only a few of which will be from the neighbourhood of x=0, where 
tan «/x“4/x.) We can avoid this difficulty by squaring each term in the 
table and working with tan?4/x (which behaves like x near x=0); 
after the interpolation has been carried out successfully for tan? (075) 
take the square root. 


Use divided differences for whole subtabulation, or alternatively to 
tabulate at x—0(-3)2-4 followed by Newton-Gregory at constant 
interval. Former method probably simpler once divided differences 
have been calculated. Check final table by differencing. Actually 
f(x)=cosh x and solution may be checked from tables, though it is 
essential that the computer should gain confidence in his own checks 
for ensuring accuracy. 


. Everett remainder same as Bessel; Stirling remainder term is 


0G? —1*) (P—2") ... P—n") pans pent) (ἢ). 


(2n-+1)! 
95 


ANSWERS AND NOTES ON EXERCISES 


Chapter Four 


1. Differentiate Bessel’s interpolation formula; use result f’= —f to 
check result by interpolation on table given. 


2. Obtained from finite difference formula by breaking down the first, 
and the first two terms, respectively, into tabular values. 

3. (ii) --ἰαὐρο 2/140; (iii) —5A3f’(O)/6 compared with —A*f’(£)/12 
for each iii of the trapezoidal rule and —h°f)(£)/90 for each 
application of Simpson’s rule. Dufton’s rule is comparable with the 
trapezoidal rule in accuracy and is remarkably easy to apply: it is use 
ful for a quick estimation of an integral from a table. 


5. Correct value is 0-86638, 
6. Correct value is 0-37558. 
7. Obtained by integration of Stirling’s interpolation formula. Second 


form for integration over a double interval obtained by adding 
formulae of first type for 


[Ὁ {{Π-{{{- {0} 


and last term obtained by suitable modification of original formula. 


= 


A h 
Ϊ Ie γ)άν τος {δ δίῳ +f tin tha th}: 
Matching is quadratic, and correction term is 
ZI fecal €:n)+ 10 χευνίξ,η) +Syvvé)}- 


Chapter Five 


1. 16—1-3x. Note that the rms residual error is about 1, and there is 
nothing to be gained from specifying the line more accurately. 


2. 5-85x**?, with rms residual error of about 0-03, 


Directory of Numerical Processes 


Polynomial approximation 

Exact matching of equidistant tabular values (see Interpolation) 
Least squares matching of discrete tabular values: §§ 5.1, 5.2. 
Continuous approximation by least squares: problem 6 (§ 5.2). 


Other approximations 


Exponential and power law approximation by least squares: 
§ 5.2. 


Interpolation 

Using divided differences based on table with variable interval: 
§ 3.2. 

Linear interpolation: § 3.1. 

Using differences of a constant interval table: §§ 3.4, 3.5: 
including throwback: ὃ 3.5. 


Differentiation 

Using differences of a constant interval table: § 4.1. 

Using tabular values only: exercise 2 (Ch. 4). 

Half-way formula using finite differences: exercise 1 (Ch. 4). 


Integration (Quadrature) 


Using ordinates at constant interval without differences: § 4.2. 
Using constant interval differences: § 4.3. | 
When the integrand has a singularity: problem 4 (§ 4.3). 
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Double and repeated integration: ὃ 4.4, exercise 7 (Ch. 4). 
Using ordinates at ‘optimum’ points: § 4.5. 
Using divided differences: exercise 4 (Ch. 4). 


Other processes 


Detection of errors in a table: example 2 (§ 2.5). 

Extraction of roots of an algebraic equation: example 8 (§ 3.6). 
Interpretation of a number: § 1.3. 

Line of ‘best fit’ for (experimental) data: § 5.1. 

Smoothing of irregular data: § 5.3. 

Summation of series: problem 5 (§ 4.3). 
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Index 
Accuracy of a number, 7 modified for throwback, 46, 
51 
Backward differences, 21 of a polynomial, 22, 23, 29 
Bessel coefficients, 42, 52 Differentiation, 56, 57 


Bessel interpolation formula, 
42, 44 
remainder, 55 
Blunders, 10, 11 


Central differences, 21 

Continuous approximation, 
90 

Continuous functions, 8 

Critical tables of coefficients, 
52 

Curve fitting, 9, 83 


D, 15 

Data smoothing, 90 

Decimal digits, 7 

ὃ, 22 

A, 20 

Detection of errors in tables, 

25, 26 

Differences, 17, 18 
backward, 21 
central, 21 
divided, 28, 29, 30 
expected fluctuations, 20 
forward, 20 
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numerical processes, 59, 80 
errors, 60 
Divided differences, 28 
Double integration, 75, 77 
Simpson’s rule, 75 
Dufton’s rule, 81 


<6, 37 
Error fan, 25 
Errors, 10 
detection, 25 
in numerical differentiation, 
57, 60 
in numerical integration, 64 
probable, 10 
random, 9 
round-off, 10, 19 
systematic, 9 
truncation, 11 
Euler-Maclaurin summation 
formula, 73 
Everett interpolation formula, 
42, 44 
Exponential approximation, 88 
Extrapolation (Newton-Greg- 
ory formula), 40 


INDEX 


Finite differences (see Dif- 
ferences) 
Floating point decimal form, 


Forward differences, 20 
Gaussian quadrature, 78, 79 


Integration, 58, 62 
divided difference formula, 


Simpson’s rule, 64, 68 
three-eighths rule, 64, 80 
trapezoidal rule, 63, 67 
Weddle’s rule, 80 
Interpolation, 8, 32, 33, 56 
Bessel formula, 42, 44 
choice of process, 44, 47 
divided differences, 34, 36 
Everett formula, 42, 44 
half-way formula, 43 
inverse, 49 
linear, 32 
Lagrange formula, 38 
Newton-Gregory formula, 
39 
Stirling formula, 43 
Inverse interpolation, 49 
Iterative process, 51 


continuous approximation, 
90 

exponential approximation, 
88 

linear approximation, 84, 
87 

polynomial approximation, 
85, 87 

power law approximation, 
88 

r.m.s. residual, 85 

Linear interpolation, 32 


μ, 42 


σ, 21 

Newton-Gregory interpolation 
formula, 39 

Newton interpolation formu- 
la, 36 

Normal equations, 86 

Notation for finite differences, 
20 

for tables, 15 
Numbers, 7 


Operational methods, 93 


Periodic functions, 16 

Pivotal values, 48 

Polynomial representation, 
3-6, 24 

Power law approximation, 88 


Lagrange interpolation formu- Proportional parts, 34 
la, 38 
Least squares, 84 Quadrature (see integration) 
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INDEX 


Random error, 9 
Round-off errors, 10, 19, 20 
Rounding-off procedure, 10 


5,15 

Semi-empirical theories, 53 

Significant figures, 7 

Simpson’s rule, 64, 68, 75 

Smoothing, 9, 90 

Stirling interpolation formula, 
43 

subtabulation, 48 
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Summation of series, 73 
Systematic error, 9 


Tables of values, 15 
Tabular interval, 15, 20 
Taylor series, 2, 53 
Throwback, 45 
Trapezoidal rule, 63, 67 
Truncation error, 3 


Weierstrass’s theorem, 4 
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Partial Derivatives P J Hilton 
| Principles of Dynamics M B Glauert 
Probability Theory A M Arthurs 
| Sequences and Series J A Green 
sets and Groups J A Green 
Sets and Numbers Ὁ Swierczkowski 
| Solid Geometry ΡΜ Cohn 
| Solutions of Laplace’s 
| Equation D R Bland 
: Vibrating Strings D R Bland 
Vibrating Systems R F Chisnell 


* Also available in cloth edition 
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LIBRARY OF MATHEMATICS 
Edited by Walter Ledermann 


The aim of this series is to provide short introductory 
text-books for the topics which are normally covered in 
the first two years of mathematics courses at Universities 
and — of Technology. Each volume is made as 
nearly self-contained as possible, with exercises and an- 
swers, and contains an amount of material that can be 
covered in about twenty lectures. Thus each student will 
be able to build ἂν a collection of text-books which is 
adapted to the syllabus he has to follow. 

The exposition is kept at an elementary level with due 
regard to modern ree. ea of aoe. Ween it is not 
feasible to give a complete treatment, because this would 
go beyond the scope of the book, the assumptions are 
fully explained and the reader is referred to appropriate 
works in the literature. 

‘The authors obviously understand the difficulties of 
undergraduates. Their treatment is more rigorous than 
what students will have been used to at school, and yet it is 
remarkably clear. 

‘All the books contain worked examples in the text and 
exercises at the ends of the chapters. They will be in- 
valuable to undergraduates. Pupils in their last year at 
school, too, will find them useful and stimulating. They 
will learn the university — to work they have 
already done, and will gain a foretaste of what awaits them 
in the future.’ — The Times Educational Supplement 

‘It will prove a valuable corpus. A great improvement 
on many works published in the past with a similar 
objective.’ — The Times Literary Supplement 

‘These are all useful little books, and topics suitable for 
similar treatment are doubtless under consideration by the 
editor of the series.’ — τ. A. A, BROADBENT, Nature 
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A complete list of books in the series appears on the 
inside back cover. 
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